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The conductivity of the two-dimensional Hubbard model is particularly relevant for high- 
temperature superconductors. Vertex corrections are expected to be important because of strongly 
momentum dependent self-energies. To attack this problem, one must also take into account 
the Mermin- Wagner theorem, the Pauli principle and crucial sum rules in order to reach non- 
perturbative regimes. Here, we use the Two-Particle Self-Consistent approach that satisfies these 
constraints. This approach is reliable from weak to intermediate coupling. A functional derivative 
approach ensures that vertex corrections are included in a way that satisfies the f sum-rule. The two 
types of vertex corrections that we find are the antiferromagnetic analogs of the Maki-Thompson and 
Aslamasov-Larkin contributions of superconducting fluctuations to the conductivity but, contrary 
to the latter, they include non-perturbative effects. The resulting analytical expressions must be 
evaluated numerically. The calculations are impossible unless a number of advanced numerical algo- 
rithms are used. These algorithms make extensive use of fast Fourier transforms, cubic splines and 
asymptotic forms. A maximum entropy approach is specially developed for analytical continuation 
of our results. These algorithms are explained in detail in appendices. The numerical results are for 
nearest neighbor hoppings. In the pseudogap regime induced by two-dimensional antiferromagnetic 
fluctuations, the effect of vertex corrections is dramatic. Without vertex corrections the resistivity 
increases as we enter the pseudogap regime. Adding vertex corrections leads to a drop in resistivity, 
as observed in some high temperature superconductors. At high temperature, the resistivity satu- 
rates at the Ioffe-Regel limit. At the quantum critical point and beyond, the resistivity displays both 
linear and quadratic temperature dependence and there is a correlation between the linear term and 
the superconducting transition temperature. A hump is observed in the mid-infrared range of the 
optical conductivity in the presence of antiferromagnetic fluctuations. 



I. INTRODUCTION 

The calculation of transport quantities in strongly cor- 
related electron systems is particularly challenging, but 
is a necessary step to make contact with a wide class of 
experiments. Even for the simplest model, namely the 
single-band Hubbard model, this is a formidable task. 
Taking up the challenge is all the more important for the 
two-dimensional case, that acts as the minimal model for 
the high temperature cuprate superconductors, 1 layered 
organic superconductors, 2 and a number of other mate- 
rials. 

Even in cases where one has a good handle on the 
single-particle Green's function, the difficulty of calcu- 
lating transport in the 2D Hubbard model stems from 
the fact that one cannot neglect the effect of vertex cor- 
rections when strong momentum-dependent correlations 
are present. Those vertex corrections are the analog of 
the self-energy, but for the two-particle response func- 
tions. When vertex corrections are not included, conser- 
vation laws can be violated and results inaccurate. In 
the case of small finite systems tractable by exact diago- 
nalization or quantum Monte Carlo calculations (QMC), 
the correlation function is directly evaluated and vertex 
corrections are not an issue. However, those results are 
more relevant for finite frequency conductivity and strong 
coupling, where correlations are mainly local. 3-8 

Consider for example the electrical conductivity for 



the two-dimensional Hubbard model. For the infinite 
system, optical and DC conductivity calculations have 
been performed without vertex corrections using Dynam- 
ical Mean-Field Theory (DMFT) 9 and Cellular-DMFT 
(CDMFT) 10 . Those calculations have also been done 
with the composite operator method 11 but vertex cor- 
rections cannot all be taken into account. For the t — J 
model, the strong coupling limit of the Hubbard model, 
a number of approaches have been used, in particular the 
extended dynamical cluster approximation 12 ' 13 , but ver- 
tex corrections 14 ' 15 have been neglected. However, recent 
optical conductivity calculations for the Hubbard model 
with the dynamical cluster approximation (DCA) took 
vertex corrections into account 16 ' 17 . The effects were 
found to be important only at high frequency. Despite 
this recent advance, the calculation of vertex corrections 
with DCA or CDMFT, considered the best available ones 
at strong coupling, 18-20 is still an open problem. 

At weak coupling, the Boltzmann equation offers a 
tractable approach that satisfies conservation laws when 
one includes scattering-in terms. For linear response, 
its variational formulation has been used for example to 
compute the effect of spin fluctuations within the self- 
consistent renormalizcd approach 21 and also to investi- 
gate the resistivity near the quantum critical point in the 
clean 22 and disordered cases. 23 ' 24 The drawback of this 
approach is that it assumes the existence of quasiparti- 
cles and that this assumption is not valid in two dimcn- 
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sions, especially near the pseudogap regime and quan- 
tum critical point. Green's function approaches that do 
not assume quasiparticles are preferable. Hence, some 
resistivity calculations without vertex corrections were 
done with the T-matrix approximation 25 and with the 
fluctuation-exchange (FLEX) approximation 26 . Other 
FLEX calculations take into account some vertex correc- 
tions due to spin and charge fluctuations 27 ' 28 . In these 
works, only the antiferromagnetic Maki-Thompson (MK) 
diagrams were included. The antiferromagnetic analogs 
of the Aslamasov-Larkin (AL) diagrams were neglected, 
claiming that the latter are negligible. We will see that in 
the DC case this is true only in the presence of particle- 
hole symmetry. A review on those calculations is given 
in Ref. 29. In Ref. 30 that takes into account super- 
conducting fluctuations, AL diagrams are taken into ac- 
count for superconducting fluctuations but they are ne- 
glected for the spin fluctuations. A recent calculation 
using field-theory methods of the conductivity at the 
quantum critical point 31 includes all vertex corrections 
but only at T = 0. The renormalized classical regime 
where a pseudogap appears has also been considered but 
focussing only on the hot spots 32 or neglecting the AL 
contribution 33 . There are also analytical results for the 
conductivity with vertex corrections using Fermi liquid 
theory 34 ' 35 . 

Despite all these results, the electrical resistivity at 
weak to intermediate coupling has not been reliably com- 
puted for all dopings and temperatures because particle- 
hole symmetry, where the AL term vanishes, does not 
generally hold and because there are regimes where Fermi 
liquid theory is no-longer valid. Fermi-liquid theory 
breaks down in the pseudogap regime, near the antifer- 
romagnetic quantum critical point and in the Ioffe-Regel 
limit. Hence, in this paper, we extend the Two-Particle 
Self-Consistent (TPSC) approach 36 " 39 to include the ef- 
fect of vertex corrections in the calculation of the resistiv- 
ity and optical conductivity of the one-band, square lat- 
tice, nearest-neighbor two-dimensional Hubbard model 
for weak to intermediate coupling. This regime corre- 
sponds to values of the interaction strength U below the 
critical value for the Mott transition. We present numer- 
ical results as examples and discuss possible links with 
experiments on cuprates. In particular, we consider the 
origin of the mid-infrared hump in the electron-doped 
materials, the Ioffe-Regel limit, insulating behavior in the 
pseudogap regime and the link between linear resistivity, 
quantum critical behavior and superconductivity. 

The TPSC approach has the following strengths that 
make it a good choice for the present purposes. In two 
dimensions, the Mermin- Wagner theorem 40 ' 41 prevents 
the occurrence of antiferromagnetic long-range order at 
finite temperature. Not many theories can handle that 
constraint. Because long-range order is prohibited, there 
is a wide range of temperatures where there are huge 
antiferromagnetic [or spin-density wave (SDW)] fluctua- 
tions in the paramagnetic state. It is in this regime that 
a fluctuation-induced pseudogap can appear. 37,42-44 The 



standard way to treat fluctuations in many-body the- 
ory, the Random Phase Approximation (RPA), leads in- 
stead to long-range order and misses this effect. The RPA 
also violates the Pauli principle in an important way. 36 
The FLEX 45,46 approximation and the self-consistent 
renormalized theory of Moriya-Lonzarich 47-49 satisfy the 
Mermin- Wagner theorem but they do not satisfy the 
Pauli principle and consistency between one and two- 
particle quantities. Strengths and weaknesses of these ap- 
proaches are discussed further in Refs. 37 and 38. Weak 
coupling renormalization group approaches become un- 
controlled when the antiferromagnetic fluctuations begin 
to diverge 50-53 . Other approaches include the effective 
spin-Hamiltonian approach 54 . The TPSC approach does 
not assume a Migdal theorem for spin fluctuations and 
Kanamori-Brucckncr renormalization of the bare inter- 
action is included without adjustable parameter. The 
conditions for the appearance of a pseudogap induced 
by antiferromagnetic fluctuations have been verified ex- 
perimentally in electron-doped cuprates. 55 In addition 
to the above theoretical considerations, the TPSC ap- 
proach has been extensively benchmarked against Quan- 
tum Monte Carlo calculations in regimes where the lat- 
ter is available. 36 ' 37 ' 39 ' 43,56-58 The agreement between the 
one-particle spectral function of the approach and QMC 
calculations in the pseudogap regime is remarkable 43 . 
The TPSC approach however fails when temperature is 
too far below that where the pseudogap appears. 

In the TPSC method, the calculation proceeds in two 
steps. At the first step, spin-spin and charge-charge cor- 
relation functions are obtained with irreducible vertices 
that are determined self-consistently. That is the ori- 
gin of the name of the approach. At the second level 
of approximation, a non-trivial self-energy that is consis- 
tent with the spin and charge fluctuations and that can 
explain fluctuation-induced momentum-dependent pseu- 
dogaps is then calculated. The charge-charge correla- 
tion function at the first level of approximation satisfies 
the /-sum rule with the Green's function at the same 
level, but it misses lifetime effects necessary to obtain 
non-trivial conductivity. What is needed is a calcula- 
tion of the current-current correlation function that in- 
cludes Green's functions dressed at the second level of 
approximation with the corresponding irreducible ver- 
tices. What has been missing up to now is an expres- 
sion for the corresponding irreducible vertices. Follow- 
ing Baym and Kadanoff, 59,60 here we use a functional 
derivative approach to obtain irreducible vertices that 
satisfy conservation laws. We check that the /-sum rule 
is then satisfied with the Green's function obtained at 
the second-level of approximation. For the conductivity, 
we show that, not only the Maki-Thompson- type con- 
tributions from spin-density wave (SDW) fluctuations, 
but also the Aslamasov-Larkin contributions arc impor- 
tant. The latter have a dramatic effect in the pseudogap 
regime. 

The paper is structured as follow. The next section 
contains the details of the methodology and is divided 
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into five subsections: II A, the model; II B, the conduc- 
tivity in linear response theory; II C, a derivation of the 
TPSC approach for the spin and charge response func- 
tions and the one-particle self-energy; II D, the conduc- 
tivity in the TPSC approach; and finally HE, a descrip- 
tion of the numerical algorithms that were used to calcu- 
late the expression given in subsection II D. Section III 
presents the results for the system considered, followed 
by a discussion and a conclusion in sections IV and V, 
respectively. Also, some useful derivations, such as those 
for the conductivity, the /-sum rule, Ward identities, are 
given in appendices along with details of algorithms that 
are of more general applicability, such as calculating re- 
sponse functions with Fast-Fourier transforms and cubic 
splines and analytical continuation of numerical data. 



II. METHODOLOGY 

We first define the model, recall the conductivity for- 
mula and introduce the TPSC method in the functional 
derivative formalism. This approach allows us, in the 
fourth subsection, to derive the conductivity formula in- 
cluding vertex corrections. The last subsection briefly 
describes the numerical algorithms that we implemented. 
Those are detailed in appendices. In this section, we use 
units where e = 1, h = 1 and a = 1, a being the lattice 
parameter. In section III, we use both those units and 
physical units to allow comparison with typical experi- 
mental values. 



A. Model 

Our model is the two-dimensional Hubbard Hamilto- 
nian in the presence of an electromagnetic field that is 
treated classically. With the usual Peierls substitution, 
we have 

H(t) -E^iv'^^^' + ^Vii. 

ija i 

(i) 

where = — rj , are the hopping matrix elements 
between the sites of the lattice, Cj a destroys a particle 
with spin a at site j and c\ a creates a particle at site i, 
A(r, t) is the vector potential, U is the on-site repulsion 
energy and n irT — c\ a Ci a is the spin a particle number 
operator at site i. Note that it is perfectly general to use 
only the vector potential to represent the field, since a 
scalar potential can always be removed using the proper 
gauge transformation. The form (1) for the Hamiltonian 
is justified by gauge invariance. Further discussion of the 
Peierls substitution may be found in Rcf. 61. 



B. Conductivity in linear response theory 

This derivation of the linear response result will allow 
us to set the notation. To obtain the conductivity, we 
first need the expression for the current operator. In the 
x direction, for example, it is given by 

where the superscript S indicates that j x is in Shrodingcr 
representation, despite its dependance on t. If we apply 
this definition to Eq. (1), keep terms up to linear order in 
the vector potential, assuming that A(r,i) = A x (r,t)x, 
with A x (r, t) varying slowly on the scale of a lattice spac- 
ing so that 

J 3 d Xij A x (r,t) « ^ {A x {v h t) + A x {vj,t)) , (3) 

where X ij IS the x component of the vector Vj — rj, we 
obtain, 

3x (r/,*) = \ J2 S * f s ( c L c l-6,<? + 4+s,* c i,°) 

Scr 

- \A x {r u t) (4<T C l-*,° + 4+5,a c l") . ( 4 ) 

Sty 

where 8 X is the projection along x of the vector 5 be- 
tween neighbors, tg is the corresponding hopping matrix 
element. For a uniform electric field we can take the vec- 
tor potential independent of position, which means that 
we need only the q = component of the current 

■Sr,\ 1 ^ e k t 

-^(^E§ c U- (5) 

ker x 

where is the dispersion relation and the number of 
lattice sites or wave vectors in the Brillouin zone. For the 
following we need to define the paramagnetic current, 

^-ArEfgcU,. (6) 

kcr 

and the diamagnetic current, 

3 d x (t) = -A x (t)±J2^ c l c ^> ( ? ) 

kcr x 

so ttet j s x {t) =3% +j d x (t). 

According to linear response theory, the frequency de- 
pendent current in response to the field is 

tix(u)) - (5'x(w)) + Xj.j. (w)A x (w) (8) 
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where j x (co) is the Fourier transform of 

j x (t) = U\t, -oo)jf (t)tf(t, -oo) , (9) 

where U(t,t') is the time evolution operator for the 
Hamiltonian H(t), Eq.(l), 

Xi.M = dte^-V (OmJUt')}) 0(t t') , 

(10) 

is the current- current correlation function and the nota- 
tion 0(t) stands for the interaction representation of the 
operator O, namely, 0(t) = e lH otQ e -^H t ^ wuere jj j s 
the Hamiltonian Eq.(l) with A = 0. Also, (...) means 
an equilibrium average, namely, in the system H . 

Since (^j x (t)j = (j x ) = 0, the equilibrium average of 
the current operator in the interaction representation is 
given by the diamagnetic term 



l(t)) = (#(*)) = (.f x (t)) 



^(*)^E^( c U, CT ). (id 

kcr 



Defining 



(k x ) 



N 



d 2 e k 



N ^ dkl 

k x 



(n k ) , 



(12) 



where n k = c^ ka Ck. a , we have 

(h(t)) = (k x ) A x (t) (13) 

or, in frequency, 

(jxM) = <**M*M , (14) 

and Eq.(8) therefore becomes 

OxM) = l(k x ) + X;w»] • (15) 

Note that, in the case where the Hamiltonian has only 
nearest-neighbor hoppings, (k x ) is proportional to the 
local kinetic energy, 8 but in general it cannot be regarded 
as such as is clear from Eq.(12). 

To find the conductivity, it suffices to relate the electric 
field to the vector potential through 



E x (ui) = + ir})A x (u) , 



(16) 



where n is positive and infinitesimal. Thus, the current 
is related to the electric field by 



t{u! + irj) 



(17) 



and finally, since (j x (uj)) — <j xx {uj)E x {u)) by definition, 
the expression for the optical conductivity in linear re- 
sponse theory is 



ct xx (uj) = 



i(uj + in) 



(18) 



In this work, we calculate only the real part of a(oj), 
the dissipative part, and the expression that we use in 
practice is 



Re a xx (uj) 



(19) 



which is derived in Appendix A from Eq.(A14) to 
Eq.(A16). If one is interested in the imaginary part, it 
can be obtained using Kramers-Kronig relations. 



C. Two-Particle Self-Consistent approach 

In the TPSC approach, one- and two-particle Green's 
functions for the Hubbard model are calculated in a non- 
perturbative way. The approach enforces conservation 
laws, key sum rules and the Pauli principle. It was 
shown, from benchmarks with quantum Monte Carlo re- 
sults 36 ' 37 ' 43 ' 56 ' 58,62 , to be accurate within a few percent 
for interaction strengths up to about U = 6t. We will 
derive the TPSC approach below, but the reader can 
also resort to Refs. 37, 38,58 and 39 for a more detailed 
discussion of the approach itself and a comparison with 
other approaches. 

In this subsection, we present the key equations for the 
theory. The following subsection contains details of the 
derivation. We use the short-hand notation 1 = (ri,ri) 
for space and imaginary time coordinates and q = (q, iq n ) 
for reciprocal space and Matsubara frequency coordi- 
nates. 

The spin and charge response functions must be com- 
puted first from the expressions, 



Xs P (q) = l f J dre iw - T <T T S*(q,r)5*(-q)) 
Xo(o) 



(20) 



fXo(<7) 



and 



Xch(q) = ^J dre*"» T (T T n(q,r)n(-q)) - (nf 



xo(q) 



(21) 



l + ¥xo(4) 



where T T is the imaginary izme-ordering operator and 
Xo(q) is the Lindhard function, given by 

Xo(q) = -2^GW(k + q)GW(k). (22) 

k 

Here, corresponds numerically to a non-interacting 
Green's function because the initial approximation for 
the self-energy, that is given below in Eq.(49), is con- 
stant and the chemical potential is adjusted to obtain 
the correct filling. This point will become clear in the 
next subsection. Note that we assume here that the sys- 
tem is paramagnetic so that the spin index will often be 
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omitted and the sum over it replaced by a factor of two, 
as in Eq.(22). In Eq.(20) and Eq.(21), the parameters 
U sp and U c h are the spin and charge irreducible vertices, 
respectively. First, U sp is defined by 



Usp %(!)> <»*(!)> 



(23) 



(this definition is used for hole doping, electron doping 
will be discussed below), so that it can be determined 
from the fluctuation-dissipation theorem 



N 



J2x sp (i) = (s z s z ) 



(n) - 2 (r^ni) 



(24) 



(n) 



U 2 



where we have used the Pauli principle (n 2 a ) = (n CT ). 
Note that all quantities on right-hand side of Eq.(24) are 
local ones. Then, once double occupancy is known, U c h 
can be determined from 



N 



E Xch{q) = (n) + 2 (n t ni) - {n)' 



(25) 



We also call Eq.(24) and Eq.(25) the local spin and local 
charge sum rules. 

Finally, with the response functions (20) and (21) we 
can obtain the one-particle self-energy, 

^{k) = Un^ 

+ J E WsvXsM + U chX ch(q)} G^(k + q) (26) 

whose form will be explained below. Expressions (20), 
(21) and (26) are the basic TPSC equations. 



1. First step: spin and charge susceptibilities 

First let us derive expressions (20) and (21) for the 
spin and charge susceptibilities. In the following, we use 
Einstein's convention for the sums (or integrals) , namely 
that an index appearing twice or more in an expression is 
summed over lattice sites and integrated over imaginary 
time. An overbar helps clarify which indices are involved. 
The approach follows the Martin-Schwinger techniques 63 
described in Kadanoff and Baym's book 60 . 

It is convenient to introduce a "source field" </> CT (l,2) 
that couples to one-particle excitations in the system. 
It allows us to easily obtain correlation functions from 
functional derivatives. The source field can be set to zero 
at the end of the calculation. The source-field dependent 
Green's function in the grand canonical ensemble is then 



given by 



Tr 



G ff (l, 2; {«£}) = - 



-^T T e-4(i)^(i,2) C ,(2) C(T (l) c t(2) 



Tr 



e -^T r e-4(i)<Mi.2>5(2) 



= - (T T <v(l)4(2)\ 



(27) 



where K = H — /j,N, fi being the chemical potential and 
N the number operator. We have also used the notation 
(. . .)^ which means that the average is taken with the 
source field turned on. Response functions can be ob- 
tained from functional derivatives of G a with respect to 
do- 1 since 



6G a {l,2;{4>}) 
5M3,4) 



We also have 



= G„,(4,3;{0})G„(1,2;{0}) 



T T ct,(3)<v(4) C t(2Ml) 



(28) 



Xsp (l,2) = (T T S z (l)S z (2)) 

= {T T K(l) - H (l)} K(2) - H (2)}) (29) 
= 2«T r n t (l)n t (2)) - (T T n r (l) % (2)» , 

where we have used spin rotation symmetry to obtain the 
last line from the previous one. For the charge response 
function, the corresponding expression is 

Xc fc(l,2) = <n(l)n(2))-<n(l)) (n(2)> 

= 2 «T r n t (l)n t (2)) + (T r n t (1)^(2))) (30) 
- Ml)) (n(2)) . 

According to Eq. (28), the last two results Eqs. (29) and 
(30) can be written as 



Xch/sp(l, 2 ) — 

_ / JG t (l,l + ;{0» JG t (l,l + ;{0}) \ 
%(2+,2) %(2+,2) ) 



(31) 



{<P}=o 



where the expressions with the plus and minus sign corre- 
spond respectively to the charge and spin response func- 
tions. Here, 1 + = {t\,t\ + e), where e is positive and 
infinitesimal. For the remainder of this section, we im- 
plicitly assume that derivatives with respect to <p are eval- 
uated at {(/)} — 0. 

To obtain integral equations for the response functions, 
one begins with 



G (J (l,l)G- 1 (l,2) = J(l-2) 



(32) 



Taking the functional derivative of this equation with 
respect to </> CT /(3, 4), taking 2—^2, multiplying on the 
right by G CT (2, 2) and summing over 2, we obtain 

^ (1 ' 2) =-G CT (l,l)^fflG CT (2,2). (33) 



MM) 



MM) 
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On the other hand, Dyson's equation in the presence of 
the field 4> reads 

G- 1 (1,2) = GS 0) - 1 (1,2)-^(1,2)-E CT (1,2), (34) 

where is the non interacting Green's function and 
X CT is the self-energy, so that 

&g-\\,2) <re g (i,2) 

= ^ ~ ^ ~ 4) ^' ~ ' (35) 

and therefore, from (33), 

<5G CT (1,2) 



MM) 



G CT (l,3)G a (4,2)<W 



Following Luttinger and Ward 64 , is a functional of 
G a and G_ CT , we find, applying the chain rule, that 



5G a {l,2) 
50a' (3, 4) 



= G < ,(l,3)G„(4,2)J (y 



This is the analog of the Bethc-Salpeter equation for 
the particle-hole channel. Defining the particle-hole ir- 
reducible vertex 



1^(1, 2; 3, 4) 



<?S g (l,2) 
^(3,4) ' 



(38) 



Eq. (37) reads 
<5G CT (1,2) 



^'(3,4) 



= G CT (1,3)G CT (4,2)5 



+ G CT (1, 1)1^(1, 2; 3, 4) ^ g( ^ } G g (2, 2) . (39) 

The spin and charge response functions in Eq. (31) are 
special cases of the more general functions 

(40) 

It is straightforward to show, from the Bcthc-Salpeter 
equation Eq. (39) and spin-rotational invariance, that 

X±(1,2;3,4) = -2G(1,3)G(4,2) 

± T cVsp (T, 2; 3, 4)G(1, 1)G(2, 2) X± (3, 4|3, 4) (41) 

where 

r ^ t (l,2) JS t (l,2) 

T cVsp (l, 2, 3, 4) - ± (42) 

and G = G t = G 4 . 



Up to now, all the results are exact. The first step 
in the TPSC method is to obtain a first approximation 
for the self-energy that will be used to obtain irreducible 
vertices. First, let us rewrite Dyson's equation Eq. (34) 
with zero source field in the form, 

G( r °)- 1 (l,3)G CT (3,2)=<5(l-2) + S CT (l,3)G <T (3,2). (43) 

Then, note that the equation of motion for the Green's 
function reads 



E 









%i + Ui 



G a (l-j,T) = S(T)5 ij 
- U{T T n ia {T)c la {r)c] a ) . (44) 



or, in compact form, 



G( r °)- 1 (l,3)G ff (3,2) = 5(l-2)-U(T T nz(l)c a (l)ct (2)) , 

(45) 

where a = —a. By comparing Eqs.(43) and (45) one 
concludes that 

S ff (l, 3)G CT (3, 2) = -U (T T n s (l)<v(l)4(2)) . (46) 

Now, comes our first approximation. We replace the 
last expression Eq. (46) by 

S CT (1, 1)G CT (1, 2) « Ug n {l)G s (1, 1+)G CT (1, 2) , (47) 

where, for the hole-doped case, 

(%(iH(i)> 



9n( l ) = 



(n t (l)) (H(l)) 



(48) 



For the electron-doped case, one replaces n a in this ex- 
pression by 1 — n a . Thus, when the lattice is bipartite, 
particle-hole symetry of the phase diagram is preserved. 
It also gives a better agreement with quantum Monte 
Carlo results when the lattice is not bipartite. Those two 
different approximation are equivalent to assume that the 
approximation (47) is made for holes when hole density 
is smaller than particle density and for particles other- 
wise. Now, substituting the definition (48) in Eq. (47), 
it is clear that one recovers the exact result for the self- 
energy Eq. (46) in the special case 2 = 1 + . The ap- 
proximation 36 ' 65 is justified if it is correct to assume that 
the four-point correlation function factorizes only when 
points 1 and 2 are different. 

Within this approximation, the self-energy can be ob- 
tained by multiplying Eq. (47) by G~ 1 (2,3) from the 
right and summing over 2, to obtain 

£W(1,2) = Ug n (l)G a (l, 1+)5(1 - 2) , (49) 

which is the self-energy ansatz that is used to define 
the Green's function G^ in the Lindhard function (22), 
from which are defined the TPSC susceptibilities (20) 
and (21). Note that once the source field c/> is turned 

off and translational invariance is restored, £^(1, 2) be- 
comes independent of position and its Fourier transform 
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is simply a constant that can be absorbed in the chemical 
potential in G^\k), so that this Green's function is in 
practice a non-interacting one. 

Given the self-energy, as well as the result 



Sgni 1 ) _ 5 mi 1 ) 



SG a (2,3) <5G_ CT (2,3) 



(50) 



valid in the paramagnetic phase, and the definition of the 
spin vertex, Eq.(42), we obtain the spin vertex 

r sp (l, 2; 3, 4) = U sp 8(l - 3)6(1+ - 4)8(1~ - 2) (51) 

where U sp — Ug^i(l). It is possible to obtain an ana- 
lytical expression for 8g^^(l) / 8G a (3,4) to compute the 
charge vertex from the definition (42) and the ansatz 
(49), but to this day the most successful approach has 
been to approximate this functional derivative as local, 
i.e. proportional to 8(1 — 3)5(1 — 4), which leads to 

T ch (l, 2; 3, 4) = U ch 8(l - 3)6(1+ - 4)5(1" - 2). (52) 

Substituting in the Bethe-Salpeter equation for the sus- 
ceptibilities, Eq.(41), we find the corresponding TPSC 
expressions for the spin (20) and charge (21) susceptibil- 
ities. They suffice to determine also U c h- 



2. Second step: improved self-energy 

Collective modes are generally less influenced by de- 
tails of the single-particle properties than the other way 
around since collective modes depend more strongly on 
general principles like conservation laws. We thus wish 
now to obtain an improved approximation for the self- 
energy that takes advantage of the fact that we have 
found accurate approximations for the low-frequency spin 
and charge fluctuations. We begin from the general def- 
inition of the self-energy Eq.(46) obtained from Dyson's 
equation. 

We start with the longitudinal channel (<fi diagonal in 
spin indices) and use the corresponding expression for 
the correlation function in terms of response function 
Eq.(28). In that case, the right-hand side of the gen- 
eral definition of the self-energy Eq. (46) may be written 
as 



E ff (l,l)G CT (I,2) = 
" 5G a (1,2) 



U 



i-.(l + ,l) 



G_ ff (1,1+) G a (1,2) 



(53) 



The last term is the Hartree-Fock contribution, which 
gives the exact result for the left-hand side in the limit 
u> — > oo. 37 The 8G a / 8<f)- a term is thus a contribution to 
lower frequencies and it comes from the spin and charge 
fluctuations. Right-multiplying the last equation by G -1 , 
replacing the lower energy part 8G a / by its general 
expression in terms of irreducible vertices, Eq.(37), and 



taking E and G on the right-hand side to be the first level 
approximations E^ and G^\ respectively, we find 



53^(1,2) = ^ (1,1+) 5 (1-2) 



UG ° (1 ' 3) ^(4,5)^(1 + ,1)' (M) 

If we expand the sum over spins on the right-hand side to 
express the irreducible vertices in terms of their spin and 
charge versions Eqs.(42) we find, after using the TPSC 
vertices, Eqs.(51,52), 

Ef)(fc) = t/n_ CT 
jj j 1 ^ 



(i) 



■U chXch (q)}G£\k + q). (55) 



There is, however, an ambiguity in obtaining the self- 
energy formula. Indeed, we can obtain an expression for 
the self-energy by using a transverse source field 4>- aa 
in Eq.(27). By taking functional derivatives of G with 
respect to this <fi, we first obtain the transverse spin cor- 
relation functions xh (1)2) = (T T S + (l)S-(2)) and x h- 

Then, using a derivation analogous to that for the longi- 
tudinal case, 38 we find 

Ef )(fc) = Un. a + H E U sp Xs P iq)G ( -l(k + q) . (56) 



During the derivation, \^ (1)2) = x h(l>2) = 

\\ S p(l,2) was used, taking spin rotational invariance 
into account. 

The two previous expressions for the self-energy clearly 
show that our approximations for the fully reducible ver- 
tex does not preserve crossing symmetry, that is the 
symmetry under the exchange of two particles or two 
holes. To improve our approximation and restore cross- 
ing symmetry, 43 we average the two expressions (55) and 
(56), which gives the final result Eq.(26) that we use in 
the rest of this paper. It turns out that this "symmetric" 
expression of the self energy gives a better agreement 
with quantum Monte Carlo results. 43 In addition, one 
verifies numerically that the exact sum rule (Appendix 
A in Rcf. 37) 



duj' 



ImE« (k,o/) = U 2 n_ a (1 



r) 



(57) 



determining the high-frequency behavior of E, is satisfied 
to a higher degree of accuracy with the symmetrized self- 
energy expression Eq.(26). 



3. Internal accuracy checks 

Our final expression for the self-energy E^ 2 \ Eq.(26), is 
in principle an improvement over the constant self-energy 
entering the calculation of the susceptibilities. But 
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clearly the approach is not one-particle self-consistent. 
All the Green's functions entering the right-hand side of 
Eq.(26) are evaluated with G^\ which has a constant 
self-energy. The advantage is that all quantities, includ- 
ing the vertices, on the right-hand side of Eq.(26) are 
computed at the same level of approximation. In fact, 
one can miss some important physics if this is not the 



case 



.37 



Apart from comparisons with quantum Monte Carlo 
(QMC) calculations, we can check the accuracy in other 
ways. For example, the /-sum rule, 

ker 

(58) 

is exactly satisfied at the first level of approximation (i.e. 
with (rikff)^ on the right-hand side) and the charge sus- 
ceptibility obtained with U c h- Suppose that on the right- 
hand side of that equation, one uses (riko-) obtained from 
instead of the Fermi function. In cases where the 
agreement with QMC calculations is good, one should 
find that the right-hand side does not change by more 
than a few percent. 

When we are in the Fermi liquid regime, another way 
to verify the accuracy of the approach is to verify if 
the Fermi surface obtained from satisfies Luttinger's 
theorem very closely. 

Finally, the consistency relation between one- and two- 
particle quantities (X and (n-j-nj,)), Eq. (46), should be 
satisfied exactly in the Hubbard model. In the special 
case where 2 = 1 + , this relation can be written as 



T 

r'A'5- N 



Tr(£G) = lim ^ £\ 



%k n T 



^ a {k)G a {k) = U{n^ ni ) . 

(59) 

In standard many-body books 66 , this expression is en- 
countered in the calculation of the free energy through a 
coupling-constant integration. In the TPSC approach, it 
is not difficult to show (Appendix B of Ref. 37) that the 
following equation 



i T r 
2 



(E( 2 >G«)=[/( W ) 



(60) 



is satisfied exactly with the self-consistent U (n-j-nj.) ob- 
tained from the sum rule (24). An internal accuracy 
check consists in ver ifying by how much ^Tr (Z^G^) 

differs from ±Tr (£< 2 )G (1) ) . A gain, in regimes where we 
have agreement with Quantum Monte Carlo calculations, 
the difference is only a few percent. 

The above relation between S and (rtj-n^) gives us an- 
other way to justify our expression for £( 2 ) , Eq. (26) . Sup- 
pose one starts from (55) to obtain a self-energy that 
contains only the longitudinal spin fluctuations and the 
charge fluctuations, as was done in the first papers on 
the TPSC approach 36 . One finds that the spin part and 
the charge part each contribute an amount U (n^n^) /2 
to the consistency relation Eq. (60) . Similarly, if we work 



only in the transverse spin channel 38,43 , we find that each 
of the two transverse spin components also contributes 
U (n^n^ /2 to |Tr (E< 2 )G (1) ). Hence, averaging the two 
expressions also preserves rotational invariance since each 
spin component contributes equally to Eq.(60). Note 
that Eq.(26) for X^ 2 ' is different from so-called Berk- 
Schrieffer type expressions 67 that contains only bare ver- 
tices and do not satisfy (Appendix E in Ref. 37) the con- 
sistency condition between one- and two-particle proper- 
ties, Eq.(59). 



D. Conductivity in the Two-Particle 
Self-Consistent approach 

To compute Re cr(u>) from the Kubo formula Eq.(19), 
we have to obtain the current- current correlation func- 
tion given by Eq. (10), which contains the q = com- 
ponent of the paramagnetic current. From the general 
expression for the current Eq.(4), we have 

n = E^) = *E te* E c l<*-',° • ( 61 ) 

( <5<x I 

Since, the actual calculation is done in Matsubara fre- 
quency instead of the real-frequency definition of Xuu 7 
Eq. (10), we use 

X jxjx (iQn) = ^J dte^-r') (T T f x {r)W)) , (62) 

where j£(r) = e TH °jPe~ TH ° and q n = 2mrT 7 with n inte- 
ger, is a bosonic Matsubara frequency. Substituting (61) 
in this expression gives 



Trf x {r)jl{r')) = - E S ^t Sl t8 2 

U5\8-2(T l<7 2 

X (Trcl(r)c^ 5l ^r)cl(T') Cl ^ a (r')) . (63) 



If we substitute the functional derivative expression for 
the correlation function, Eq.(28), in this equation, we 
obtain 

(T T f x (n)f x (r 2 )) = E ***** %l'n,v] 

(64) 

where have used the notation 1 = (ri,Ti) and V = (ri — 
5\ , n) . Note that the first term on the right-hand side of 
Eq. (28) does not contribute to the sum since = 0. 
Now, summing over spin indices and using spin rotational 
invariance, we obtain 



Trf x {Tl)f x {T2)) = - E <^ 2 *+( 2 '' 2 ' ' 

rir 2 <5i(52 

(65) 

where we have used the definition Eq. (40) for the gen- 
eralized susceptibility \ + . 
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The most general way of thinking about the last result 
is that it comes from a functional derivative of the current 
with respect to an applied vector potential representing 
the electric field. The remaining of the derivation, that 
we give in the rest of this subsection, is based on the 
idea that we should evaluate this functional derivative in 
a systematic way to obtain a conserving approximation, 
as shown by Baym 59 . We will keep working with the 
source field 4> a but it is useful to remember that within 
simple prefactors, it is equivalent to working with the 
vector potential as the source field. 

The susceptibility \+ is defined in terms of a functional 
derivative in Eq. (40). That functional derivative leads 
to a Bethe-Salpeter like equation (36) that contains two 
different types of terms 

X+ (1,2|3,4) = -2G(1,3)G(4,2) 

The product GG comes from the explicit dependence of 
G -1 on the source field while the last one comes from 
the implicit dependence of the self-energy on that source 
field. The product GG is what leads to the "bubble" 
in standard conductivity calculations. The other term is 
the vertex correction. 

Up to now, everything is exact in the present subsec- 
tion. To work within the TPSC approach, it suffices to 
use everywhere the results obtained at the second level 
of approximation, namely G& and H„ since at the first 
level of approximation the conductivity is infinite. Our 
explicit formula for Eq. (26) is a functional of G^jP . 
We can thus write 

6X?\l,2) 5Yg\l,2) 5G { P {1,1) 



MM) ^(3,4) <5<M3,4) 



(67) 



Expanding the sum over spin, using the chain rule, spin 
rotation invariance and the definition of X+ Eq. (40) we 
obtain 



6^(1,2) 6Y%> (1,2) 
<5<£ CT (3,4) <ty-„(3,4) 



(2), 



1 (S^\l,2) 



* \6G { P(3,4) <5G^(3,4) 

where x+^ is computed with G^\ Substituting in the 
exact expression for X+ Eq. (66), we find 

X f (1,2|3,4) = -2G< 2 >(1,3)G< 2 >(4,2) 

\6G { P(3,4) 5GV(3,4)J 

xxW(3,4|3,4). (69) 

If Si 2) had been a functional of G^ we would have 
had an infinite series to sum. Instead, here the series 



<raS 2) (l,2)' 



xi 1} (3,4|3,4) (68) 



ends as we now show. First note that, from Eq.(41) and 
the approximation (52) for the charge vertex, 

xV 1) (3,4|3,4)--2GW(3,3)G«(4,4) 

+ G« (3, 5)G« (6, 4)r« (5, 6, 7, 8) X « (7, 8|3, 4) 

= -2G (1) (3,3)G (1) (4,4) 

+ U ch G^(3, 5)G«(5, 4)xV 1) (5, 5+|3, 4) . 

(70) 

From this expression, only the first term has to be in- 
cluded in our calculation because the second term will 
not contribute to the q = current- current correlation 
function. This is because is local in space and local 
vertex corrections do not contribute to the uniform con- 
ductivity. To see why, one can represent the correlation 
function Eq. (65) as a series of diagrams with current ver- 
tices at their ends and vertex functions inserted between 
them. Whenever a vertex function in a diagram does not 
depend on q, a bubble is closed with the local vertex on 
one end and the current vertex on the other. Since a cur- 
rent vertex is an odd function in space and the product 
of the Green's function is even because the bubble does 
not carry any momentum, the integral vanishes. 

Inserting the first term of x+\ Eq. (70), into x+\ 
Eq. (69), we obtain 



xf (1,2|3,4) 



-2G (2) (1,3)G (2) (4,2) 

(2)/T <^ r^(2). 



x G (1) (3,3)G (1) (4,4). (71) 

The last step is in principle straightforward, but very 
tedious. We must obtain an expression for the functional 
derivatives in parenthesis in this expression. The self- 
energy Eq. (26) in real space is 

S^(l,2) = C/G«(l,l+)<5(l-2) 

+ % {3U spXsp (2, 1) + U chXch (2, 1)] G«(l,2) . (72) 



so that 

5y£ ) (1,2) 
<5G«(3,4) 

+ % [3U spXsp (2, 1) + U chXch (2, 1)] 6(1 - 3)6(2 - A)6 a 



US(1 - 2)6(1 - 3)6(1+ - 4)<5_ ff , CT , 



U 



+ ^GW(1,2) 



3U, 



Sp 



S Xsp (2,l) 
6G$(3A) 



+ U ch 



SXch(2, 1) 
6G^(3,4)\ 



(73) 



In this expression, the terms involving the functional 
derivatives of T sp and T c h have been omitted because 
they do not contribute to the conductivity. Fundamen- 
tally, what is needed is sr °^ /ch } where A x is the vector 
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potential. Since r sp and T c h are local, this is the correla- 
tion function of a local operator, that is thus even under 
parity, and the current operator, that is odd. The q = 
component of this correlation function that enters the 
current- current correlation function, Eq.(62), therefore 
vanishes. 

Now, we need an explicit expression for 
6Xch/sp(2, 1)/ (5G'£ 1 ,' ) (3, 4) to know the vertex correc- 
tion Eq. (73) completely. Let us start with Xch- One can 
obtain an expression for this function by taking 2 = 1+ 
and 3 = 4+ in the expressions for x+/- Eq. (41). But 

I 



this expression has been simplified using spin rotational 
invariance. If we separate the spin contributions, we 
have 

Xch (2, 1) - -G„(2, l+)G a (1, 2+) - G_ ff (2, l+)G_ ff (l, 2+) 
+ l -U ch G a (2X)G a {%2+)xch{%l) 
+ l -U ch G- a % 2+)G_ CT (2, 2+) X ch{% 1) 

(74) 

so that 



SXch(2,l) 
5G a ,(3,A) 



= -5(2 - 3)5(1+ - 4)G^(1, 2+) - G a ,{2, 1+)5(1 - 3)5(2+ - 4) 

+ \u ch 5(2 - 3)G a ,(A, 2) X ch(4, 1) + l -U ch 5{2 - A)G a ,(2, 3) Xc/l (3, 1) 



^ S Xch (2,l) 1 
JG„/(3,4) 2 



+ -U ch G a (2, 2+)G„(% 2+) "JTTo^ + ^ ch G_ a (2, 2+)G_ CT (2, 2+) , (75) 



^Xch(2, 1) 
<5G CT ,(3,4) 



where, again, the functional derivatives of U c h have been omitted, for the same reason as in (73). Using again spin 
rotational invariance, this gives 



SXch(2,l) 
SG a ,(3,4) 



-(5(2 - 3)5(1+ - 4)G(1, 2+) - G(2, 1+)5(1 - 3)5(2+ - 4) 
+ \u ch 5(2 - 3)G(4, 2) Xcft (4, 1) + \u ch 5(2 - 4)G(2, 3) Xc h(3, 1) + tf ch G(2, 2+)G(2, 2+) ^|^y • (76) 



By Fourier transforming this equation with respect to 2, we obtain an algebraic equation that is trivial to solve. 
Fourier transforming back the result, we obtain 



5G*'(3,4) l 



-(2, 3) G(4, 3) [ - 5(1 - 4) + \u ch Xch (4, 1) 



L Xo 



+ : , l h (2,4) G(4, 3) \- 5(1-3) + \u ch Xcfe (3, 1) . (77) 
1 + ^Xo L 2 



where 



1 + Yxo 1 ' J ^ 1 + % 



+ ^tXo(q) 



(78) 



The expression for 5xsp(2, l)/5G a i(3, A) is obtained by simply replacing U c h by —U sp and Xc/i by Xsp in the right-hand 
side of expression (77), 



<W2, 1) 



1 

T^ 



(2,3)G(4,3) 



<5(l-4)--[/ spXsp (4,l) 



+ 



1 

"Us 



■-Xo 



-(2,4)G(4,3) 



5(l-3)--U spXsp (S,l) 



(79) 
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Substituting Eqs. (77) and (79) in our expression for the vertex, Eq. (73), we obtain 

""' ' ! ' 2) ~ US(1 - 2)5(1 - 3)5(1+ - 4)5_ ctV + J [3U spXsp (2, 1) + U chX ch(2, 1)] 5(1 - 3)5(2 - A)S aa > 



5G«(3,4) 



U 



G (1) (1,2)G (1) (4,3) 



3U S 



I _ (2, 3) [5(1 - 4) + \u ap Xflp (4, 1)] + (2, 4) [5(1 - 3) + \u sp Xsp (3, 1)] 



+ U ch ( ^ (2,3) [5(l-4)-^ efeXeft (4,l)]+ ^ (2,4) [5(1 - 3) - ^ Xch (3, 1)] 

\ ~r" 2 XO 2 XO 

and therefore, 
5£<- 2) (l,2) 5ES 2) (1,2) 



(80) 



5G«(3,4) 5G«(3,4) 



-^G (1) (1,2)G (1) (4,3) 



3C/ S 



(75(1 - 2)5(1 - 3)5(1+ - 4) + ^ [3t/ spXsp (2, 1) + U chXch (2, 1)] 5(1 - 3)5(2 - 4) 

O 

p ( _ L & 3) _ 3) + ^ Xsp(3, 1} ] + 1 _ L & 3) _ 5) + ^ Xsp(3 ' T) ] ) 

V ' 2 XO * 2 XO "" / 



+ tf<* L,^ (2,3) [s(l-4)-^U chXch (i,l)] + * (2,4) [5(1 - 3) - \u ch Xch (3, 1)] ) 

\ ^ "T" 2 XO - ^2 XO ~ ~ J 



(81) 



We now have everything to write down an explicit form for the susceptibility X + m (71), 
xV 2) (l,2|3,4) = -2G< 2 )(1,3)G (2) (4,2) -2C/G (2) (1,1)G (2) (1,2)G (1) (1,3)G (1) (4,1+) 

-^G( 2 >(1,1)G( 2 >(2,2) [3C/ spXsp (2,T) + C/ c ^(2,T)]GW(l,3)G«(4,2) 



U 



G (2) (1,1)G (2) (2,2)G (1) (1,2)G (1) (4,3)G (1) (3,3)G (1) (4,4) 



3U, 



sp 



1— (2,3) \- 6(1-4) -lu spXsp (i,l) 



-Xo 



+ 



-Xo 



-(2,4) -6(1-3)- -U spXsp (3,l) 



+ U ch 



1 



1. 



(2,3) -J(l-4) + -C^,Xcfc(4,l) 



1 



1 



-(2,4) -5(l-3) + -t/cfcXch(3,l) 



(82) 



1 + ^Xo 

To evaluate the current-current correlation function entering the conductivity, it is clearly necessary to go to Fourier 
space. Inserting then this last result in Eq. (65) for the current-current correlation function and using 

N 



and 



1 + 



U, 



sp 



Xsp(q) 



l 



u. 



2 xo(q) 



dk x 

1 U <* < \ 1 
l-^Xc,( 9 )= 1 + ¥xoW , 



(83) 



(84) 



we finally obtain for the Fourier-Matsubara transformed expression at q = 0, 

2 



Xj^A^Qn) = 



U (T 



-2T 



k 



U (T 



4 ( n I S G {2 \k 1 )G^(h + iq n )G^(k 2 )G^(k 2 + ^n)^(h)^-(k 2 ) [3U spXsp (k 2 - h) + U chXch (k 2 - h)} 



fei ko 



2 VAT 



dc 



<9e k 



^ W {kl) W (fc2)G(2) (fcl)G(2) +^) G(1) ^ G(1) ( fc 2 +*9n) (fc 2 + qi + iq n ) + G« (k 2 - qi ) 



X G {l) (k 1 +q l +iq n )\3U s 



1 1 ! I 



1 - % Xo(9i) 1 - ^Xofe + iq n ) ' ' 1 + ^Xo(tfi) 1 + Y*o(<Zi + «<Zn) 



(85) 



In this expression, we use the compact notation k + iq n = (k, ik m + iq n ). Note that the second term on the right- 
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hand side of Eq. (82) does not contribute to Xj x jx because 
it has a local vertex. 




FIG. 1. Schematic representation of the various terms in 
the current- current correlation function Eq.(85). Simple 
fermion lines represent the Green's function G^ 1 ' (in practice 
it has the same form as the bare Green's function), double 
fermion lines are dressed Green's function G^ 2 \ the simple 
wiggly line in the last diagram on the first line is the func- 
tion [3U Bp Xs P + UchXch], the double wiggly lines in the third 
and fourth diagrams on the right-hand side are the function 
1/(1 — UspXo/2) and the zigzag lines in the last two diagrams 
are the function 1/(1 + [7 c hXo/2). Those diagrams serve only 
as a graphical representation of Eq.(85), they are not obtained 
in a perturbative way. 

Eq. (85) is our final result for the current-current corre- 
lation function in Matsubara space, including both bub- 
ble and vertex corrections. It is useful for the intu- 
ition and to verify the Fourier transforms to represent 
it schematically, as we have done in Fig. 1. We stress 
that this representation is not the result of a perturba- 
tive calculation. It exists merely because we are working 
with Green's functions and response functions. This rep- 
resentation helps noticing that an analogy can be made 
between the different contributions to Xj^j* an d the dia- 
grams considered in the theory of paraconductivity. 68 In 
the latter theory, one considers the effect of the super- 
conducting fluctuations on the conductivity of the normal 
state, while in our case, the bosons exchanged are spin 
and charge fluctuations. The diagram with a single bo- 
son propagator on the first line of Fig.l is the analog of 
the Maki-Thompson (MT) contribution, 69-72 while the 
diagrams on the second and third lines are analogs of 
the Aslamasov-Larkin (AL) contributions. 73 However, it 
is important to note that, because our approach is not 
perturbative, some electron propagators are at the sec- 
ond level of approximation and some are at the first level, 
and the boson propagators are the susceptibilities com- 
puted with the renormalizcd spin and charge irreducible 
vertices. 

In the system considered in this paper, the relevant 
collectives modes are magnetic fluctuations with a wave 



vector close to Q = (tt, tt). First, when those fluctuations 
become strong, they scatter quasiparticles, which has a 
dramatic effect on the single particle spectrum and there- 
fore the conductivity obtained from the bubble alone. 
But the magnetic fluctuations also lead to important ver- 
tex corrections because they correlate the regions of the 
Fermi surface that are connected by wave vectors close to 
Q. For the MT diagram in Fig. 1, in the DC limit, the 
exchange of such a fluctuation between the particle and 
the hole created by the field causes the pair to be scat- 
tered from a region of the Fermi surface where it carries 
a positive current along the direction of the field, to a re- 
gion where it carries a negative current. This correlation 
between currents flowing in oposite directions degrades 
the DC conductivity. In the case of the AL diagrams in 
Fig. 1, a particle-hole pair creates another one via two 
magnetic fluctuations. If they have a large correlation 
length, this will tend to correlate currents on large dis- 
tances. However, the particle-hole pair created by the 
fluctuations can carry a current that is either positive or 
negative. Keeping in mind that, when quasiparticles cre- 
ate or absorb magnetic fluctuations, their velocity along 
the direction of the field changes sign, one finds that the 
first diagram in the second line of Fig. 1 correlates cur- 
rents flowing in the same direction, while the second dia- 
gram correlates currents in opposite directions. In addi- 
tion, the former contribution is large if the single particle 
spectral density is large below the Fermi level in regions 
connected by wave vectors around Q, while the latter is 
large if the spectral density is large above the Fermi level. 
When both processes are summed up, there will there- 
fore be a net effect on the conductivity only if there is an 
asymmetry in the relevant parts of the density of states. 



E. Calculation algorithms 

Let us recall that we have used the short-hand no- 
tation k = (k, ik n ), with k = (k X7 k y ), in the general 
expression for the current- current correlation function 
Eq. (85). Thus, each sum over a tri- vector k is over the 
two-dimensional Brillouin zone and over Matsubara fre- 
quencies. The second and third sums are therefore six- 
and nine-dimensional, respectively. The six-dimensional 
sum would be extremely long to do in a direct way for 
relevant system sizes and temperatures. For example, 
for a finite system of 512 x 512 sites with about 4000 
frequencies, which would allow to go down to temper- 
atures about T = O.Oli (without finite size effects), it 
would take of the order of thirty years to compute one 
frequency iq n , if 10 9 terms are summed per second. Of 
course, this is using pure brute force, when all wave vec- 
tors and frequencies are summed, which is not necessary 
in practice. In the case of the nine-dimensional sum, it 
would take about 40 billion years to compute a single 
frequency in the same conditions. If one was to keep the 
direct summation approach, but optimizing the proce- 
dure by using a very efficient adaptive scheme, it would 
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still be extremely hard to do the six-dimensional sum 
to obtain, for example, 100 frequencies in a reasonnable 
time, namely, of the order of a few days. In the case of 
the nine-dimensional sum, it is obviously impossible to 
do this way. 

One thus has to resort to a numerical approach com- 
pletely different from direct summation to succeed in cal- 
culating all the terms of Eq.(85) for a useful number of 
frequencies iq n . The main tool that we use to make this 
calculation possible is the fast Fourier transform (FFT) , 
which changes the scaling of Fourier transforms (FTs) 
from iV 2 to NlogN. But other mathematical and nu- 
merical tricks are also necessary to make the calculation 
both fast and precise. Precision is critical here since we 
have to numerically perform the analytical continuation 
of the computed Matsubara current-current correlation 
function and numerical analytical continuation is inher- 
ently an ill-conditionned problem. For this analytical 
continuation procedure, that produces the real frequency 
optical conductivity from the Matsubara current-current 
correlation function, a maximum entropy algorithm is 
also developpcd to maximize the accuracy of the result. 
The calculation algorithms for Xj^jxi^Qn) and our analyt- 
ical continuation algorithm are summarized in the next 
two subsections and detailed in appendices C to F. 



1. Fast Fourier transforms, cubic splines and asymptotic 
expansions 

The key property that allows us to compute Eq.(85) 
is that some of its sums are convolutions. Since the con- 
volution of two functions can be written as the Fourier 
transform of the product of the of their Fourier trans- 
forms, we can use fast Fourier transforms (FFT) to do 
those convolutions in a very efficient way. However, FFTs 
are discrete transform, while some of the transforms we 
have to do are continuous ones. In the case of the spatial 
Fourier transforms, we can use a finite system with peri- 
odic boundary conditions, so that all the transforms are 
discrete and FFTs can be used directly. Since we work 
at finite temperature in two dimensions, all correlation 
lengths are finite and we can use a system large enough 
to reach interesting regimes. The system size we use in 
this work is 512 x 512 and the lowest temperature we 
reach is 0.0082;, which corresponds to 32AT \it~ 0.35eU. 
At this temperature, the thermal Dc Broglie wavelength 
is about 100 in lattice units, so that no finite size effect 
are seen unless the magnetic correlation length becomes 
large. 

In the case of imaginary time, if it is discretized, then 
the Fourier transform will be periodic in Matsubara fre- 
quencies, which is unphysical. In fact, the lowest frequen- 
cies will be acceptable, though not very precise, but the 
precision will decrease rapidly with increasing frequency, 
so that the high frequencies will be completely wrong. 
To overcome this problem, we use a cubic spline to inter- 
polate the function between the discrete imaginary time 



points and we do a continuous Fourier transform on the 
spline. This technique has already been used in the con- 
text of dynamical mean field theory calculations. 74 How- 
ever, to our knowledge, it has not been pointed out that, 
because the spline is only twice continuously diffcrcn- 
tiable, after integrating by parts three times the Fourier 
integral, we are left with an expression containing a dis- 
crete Fourier transform (DFT) that can be done with 
a FFT. Using the more general formula derived in ap- 
pendix E, the final formula for the Fourier transform of 
the cubic spline of an imaginary time function g(r) is 

g(icj n ) = / drg(r)e^ 
Jo 

N T] 

drSjiry^ 

j = l J Tj-l 

= ff(0) - e^g((3) + S[(0) - S' N {[3) 
iuj n (iu n ) 2 

(iuj n ) 3 

i „iw„Ar N-l 

+ 1^1 V S (3) e 4 ^ 

(86) 

where A" is the number of intervals in the discretized 
imaginary time between and j3, At is the size of an 
interval, Sj(t) is the cubic polynomial in the j interval, 
Sj(t), S"(t) and are, respectively, the first, second 
and the third derivative of Sj(r). If g is a fcrmionic 
function, e lLOnl3 = — 1, while e JLJ " /3 = 1 if it is bosonic. 

When using formula (86), we work explicitly with the 
high frequency expansion up to the 1/ (iui„) 3 term, which 
can be very useful, as will be explained shortly. Up to 
which term this expansion will be exact depends on how 
the spline is defined. As shown in appendix E, there are 
two conditions defining the spline that can be chosen de- 
pending on the information available. If the derivatives 
at the boundaries g'(0) and g'(f3) are known, one can fix 
S[(0) and S' N (f3) to those values. That is what we do in 
our calculations. Otherwise, one has to use exact results 
for the second and third high frequency coefficients of g, 
i.e., the second and third moments of its spectral func- 
tion. Another important point about Eq.(86) is that the 
last term containing the DFT only contributes to the low 
frequencies because of its factor l/(«ix>„) 4 , therefore the 
error at high frequency that comes from the discretiza- 
tion of r rapidly vanishes. 

Given the above considerations, it should now be clear 
that we can compute convolutions both in wave vector 
and Matsubara frequency using fast Fourier transforms. 
Before the calculation of Xj x jx> Eq.(85), this technique 
is used to compute the Lindhard function Eq.(22) and 
then, the self-energy Eq.(26). The computation of those 
functions with FFT is relatively straightforward to imple- 
ment, although some care must be taken in the definition 
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of the splines. Those details are given in appendices C 1 
and C2. In the case of Eq.(85), we only seek the q = 
component of Xj x j x (<1> *<7n) so that, in each term, the last 
sum over the Brillouin zone to be performed is a simple 
sum and FFTs are not useful there. Regarding the other 
sums it is not obvious which one can be put in the form of 
a convolution, except for the bubble which has the same 
form as the Lindhard function. In the case of the sec- 
ond term of Eq.(85), represented by the second diagram 
(Maki-Thompson) on the right-hand side of Fig.l, it is 
possible to put all those sums in the form of convolution, 
and thus to compute all of them with FFTs. This means 
that, all the external frequencies iq n can be obtained at 
once, when the last FFT is performed. The scaling of this 
calculation with system size and inverse temperature is 
thus the same as for the bubble. That is not the case 
however of the last term of Eq.(85), represented by the 
last four (Aslamasov-Larkin) diagrams in Fig.l. For this 
term, each frequency has to be calculated separately and 
it is the calculation time of only one frequency that has 
the same scaling as the bubble. Calculating this term for 
100 values of iq n therefore takes hundreds of times the 
calculation time of all frequencies of the bubble. The cal- 
culation of the bubble is nevertheless very quick, of the 
order of a few minutes at most for our 512 x 512 system 
at low temperature. 

Note however that a quite large amount of analytical 
work is needed to put the last two terms of Eq.(85) in 
a form suitable for computation with FFTs. This work 
involves some transformations and a certain number of 
sums over Matsubara frequencies must be performed ex- 
actly. The details of those transformations and analytical 
calculations are given in appendix C 3. 

Still, it is impossible to calculate the third term of 
Eq.(85) for thousands of frequencies, namely the num- 
ber of frequencies below the cutoff used in the calcula- 
tion. Therefore, to reduce the number of frequencies to 
be calculated, we have used a non-uniform Matsubara 
frequency grid in which the frequency spacing increases 
with frequency magnitude. We give the definition of this 
grid in appendix D. We have verified that the function 
Xj x j x evaluated on this grid has enough information 
for the analytically continued conductivity Re <j(lj) to be 
converged. 

Using the formula (86) for our functions in Matsub- 
ara frequencies has the advantage that we always have 
at our disposal their high frequency expansion. Thus if 
we want the inverse Fourier transform (from Matsubara 
frequencies to imaginary time) of a more complex func- 
tion, for example the product of two functions of the form 
(86), we can always use the asymptotic form to perform 
the sum over Matsubara frequencies up to infinity using 
contour integrals in the complex plane and the residue 
theorem. In a lot of cases, this is necessary since the 
sum does not converge otherwise. As will be explained 
in the next subsection, precision is important in our cal- 
culations since we have to extract Re <j(u>) from the Mat- 
subara function Xjxjx and analytical continuation is 



a very ill-conditioned problem. 

Finally, as one will notice from appendix C 3, the calcu- 
lation of Eq.(85), especially the second and third terms, 
contains many steps. The debugging part of the work is 
therefore considerable. To make sure the formula (85) is 
implemented correctly, we have compared its brute force 
implementation, term by term, with its fast one given in 
appendix C3. Here, brute force means that it is coded 
exactly as written in (85), with one exception that will 
be explained shortly. Therefore, we compare the results 
of two calculations that are completely different numer- 
ically, but mathematically equivalent. Of course, those 
verifications can only be done for very small system at 
very high temperature, since otherwise the brute force 
calculation is impossible to do in a reasonable time. But 
even then, the third term of (85) takes too much time to 
compute with brute force. In that case, the procedure is 
to first verify that the sum over fc 2 for some random val- 
ues of q n and q\ gives the same result with the two imple- 
mentations and then to replace this sum in the brute force 
version by its fast implementation. After that step we 
compare this modified brute force implementation with 
the fast implementation of the whole third term. 



2. Analytical continuation 

Expression (85) gives the current- current correlation 
function in Matsubara frequencies while we need it in 
real frequency to compute the conductivity from Eq.(18) 
or Eq.(19). Thus we need a reliable analytical contin- 
uation method. To do so, we use a maximum entropy 
approach. This kind of analytical continuation proce- 
dure is often used to extract real frequency results from 
imaginary time quantum Monte Carlo data using some 
information known a priori such as sum rules and a de- 
fault model which contains some known properties of the 
expected function 75 . This information is included in the 
algorithm in the form of constraints or in the entropy def- 
inition. This kind of approach is well suited for quantum 
Monte Carlo since the results are usually in discretized 
imaginary time and the amount of noise can be impor- 
tant, so that any additional information is welcome. In 
our case, the original data is in Matsubara frequencies 
and the noise level is much lower since it comes only from 
finite precision rounding errors accumulating throughout 
the calculation. However, our calculation of Xj x j x con- 
tains many steps so that the final result may have an 
amount of noise such that Pade approximants will not 
work. Those approximants are very sensitive to noise 76 
and their reliability is not very good except at very low 
temperature. In any case, we have noticed that they give 
unstable results as a function of temperature for our data 
while our maximum entropy procedure is in general sta- 
ble. 

We give here the main features of our approach. The 
details are given in appendix F. As usual we minimize the 
function \ 2 — a S, where S is the entropy, a, a weighting 
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parameter for S while \ 2 is the quadratic error between 
Xj^jxii-Qn) and the quantity Xjj{iq n ) computed from the 
spectral representation of Xj x j x (*<Zn) with a trial real fre- 
quency conductivity Re a(co). First, an accurate nu- 
merical integration scheme has to be chosen to compute 
Xjj(iq n ) from Re <j(w) known on a fixed discrete grid in 
real frequency u. We use a cubic spline to approximate 
Re a(ui). Since we work in Matsubara frequencies in- 
stead of imaginary time, the spectral form is simple and 
can be integrated analytically if Re <j(lu) is approximated 
by a piecewise cubic polynomial function. Also, because 
we want to integrate Re <j(ui) in the interval [0, oo, we 
integrate the low frequency part with respect to uj and 
the high frequency part with respect to 1/u). To do so, 
we use a spline cubic in uj for low frequencies and cubic 
in for high frequencies. As for the choice of grid, 
we take it to be uniform in u> for the low frequency part 
and uniform in l/u> for high frequencies. This choice of 
grid ensures that the matrix for the spline linear sys- 
tem is well conditioned, and keeps the number of values 
Re <j(ujj) reasonable, so that the minimization procedure 
is not too heavy. The integration using the spline turns 
out to be very accurate compared to a simple piecewise 
linear approximation. In our tests with a well defined an- 
alytical form for Re <t(uj) for which Xjxjxi^ln) could be 
computed very accurately with an adaptive integration 
routine, the relative precision was typically five orders 
of magnitude smaller with the spline than the piecewise 
linear approximation. The reason why we have to use a 
very accurate integration method is that the relative dif- 
ference between Xjj{iq n ) and Xjxjxi^ln) has to be very 
small, typically < 1CT 5 in our calculation, for Re er(cj) 
to be converged in the optimization procedure. The pre- 
cision of the numerical integral for a given Re <j(u>) has 
clearly to be smaller than this required relative error on 
Xjj(iq n ) for the result to be reliable. 

The fact that the spline is integrated analytically has 
the great advantage that low temperatures are not more 
difficult to handle, while it is the case with standard 
numerical integration because the integrand becomes 
sharper as temperature decreases. Note that all those 
precision issues are important if one is interested in quan- 
titative results. For instance, we want to obtain the resis- 
tivity as a function of temperature, but since the results 
at different temperatures are numerically completely in- 
dependent, the quantitative aspect becomes crucial. If 
one is only interested in the shape of Re <j(uj) at a given 
temperature, simpler and cruder approximations can be 
sufficient. 



III. NUMERICAL RESULTS 

Before we show the results of this work, let us recall 
some important former TPSC results. First, the theory 
respects the Mermin- Wagner theorem, so that no phase 
transition occur at finite temperature. However, with 
proper values of U and hopping parameters, antiferro- 



magnetic correlations are present up to very high tem- 
peratures around half-filling. For example, with U = 6t 
and nearest neighbor hopping only, for dopings smaller 
than p c = 0.205 a crossover to a renormalized classical 
regime appears. This regime appears when ksT » huj sp , 
where oj sp is the characteristic frequency of the antifcr- 
romagnetic fluctuations, i.e., the frequency at which the 
imaginary part of the spin correlation function x" P ( w ) 
is maximum. In this regime, the antiferromagnetic cor- 
relation length has the form £ sp oc exp(C/T), where 
C has a very weak temperature dependence. There- 
fore, at a certain temperature T* , £ sp becomes larger 
than the single-particle thermal De Broglie wavelength 
£th = tivp/ '(irksT). When this happens, the parts of 
the Fermi surface that are connected by the antiferro- 
magnetic wave vector, called the hot spots, are strongly 
scattered by the magnetic fluctuations and eventually de- 
stroyed, producing a gap in those regions of the Brillouin 
zone. 37 ' 62 However, before the correlation length becomes 
infinite, which is the case only at T = in two dimen- 
sions, there is still spectral weight at the Fermi level and 
thus no real gap exits, but what is observed instead is a 
pscudogap, namely, a depression of the density of states 
at the Fermi level. Therefore, the crossover temperature 
to the renormalized classical regime T* can also be called 
a pseudogap temperature. When T = 0, long-range spin- 
density wave (SDW) order exists for p < p c and thus p c 
is an quantum critical point (QCP). Depending on band 
parameters and doping, this SDW state can be commen- 
surate or incommensurate. Usually, it is commensurate 
close to half-filling and a transition to incommensurate 
appears at a certain doping. 77 

Benchmarks of TPSC calculations were made against 
quantum Monte Carlo (QMC) results for quantities such 
as the spin- an charge structure factors, 36 the quasi- 
particle renormalization factor and the imaginary time 
Green's function, 62 the finite frequency spin susceptibil- 
ity, the double occupancy and the one-particle occupa- 
tion number, 37 and finally, the one-particle the spectral 
weight 43 . Those benchmarks were made in the weak to 
intermediate coupling regime for a large range of doping 
around half-filling and for temperatures where no finite 
size effect are seen in the QMC results and, at finite dop- 
ing, when the sign problem is not too strong. In general, 
a good quantitative agreement is obtained for all quan- 
tities at a coupling U = At. For some quantities such as 
the spin structure factor, the agreement is almost perfect 
above T* for couplings up to U = 8t. Since the TPSC 
approach has some mean field aspects 78 coming from the 
ansatz Eq.(49), it slightly overestimates T* (see Fig. 7 
of ref. 37), but the qualitative behavior just below T*, 
when the spin correlation length grows exponentially, is 
still very well reproduced. The spectral weight at half- 
filing in this regime is also very well reproduced by TPSC 
calculations. 43 This is the regime where precursors of the 
antiferromagnetic bands are formed and a pseudogap ap- 
pears in the spectral weight. 

The results we present in this section are for the 
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FIG. 2. Contributions to the zero Matsubara frequency value 
of the current- current correlation function compared to the 
sum-rule value — (k x ) for p = 0.17. 



equation is satisfied to a relative accuracy of 10~ 7 . By 
increasing the cutoff, the accuracy can be increased at 
will. We have reached an accuracy of 10~ 10 . The sep- 
arate contributions of the different terms of expression 
(85), also represented schematically in Fig.l, are shown 
in Fig. 2 as a function of temperature for 17% doping, i.e. 
on the left side of the quantum critical point. As one 
would expect, the bubble contribution is dominant at 
high temperature, although the first vertex correction is 
not negligible. At low temperature, in the renormalized 
classical regime, the two terms give comparable contri- 
butions. The contribution of the third term in Eq.(85), 
the Aslamasov-Larkin-likc diagrams in Fig.l, vanishes at 
all temperatures. As will be seen in the next subsec- 
tions, despite this vanishing contribution of this term to 
XjxjxiiQn — 0), i.e. to the integral of the real part of the 
conductivity, this term contributes in a non-trivial way 
both to the DC and the finite frequency conductivity. 



one-band Hubbard model with nearest-neighbor hopping 
only. All numerical examples are for U = 6i and various 
dopings and temperature. We begin by showing the ac- 
curacy with which the /-sum rule is satisfied. We then 
give a few typical examples of the frequency dependent 
conductivity. The last subsection will focus on the tem- 
perature dependent resistivity for various dopings. 

A. f-sum-rule 

Although our expression for the conductivity Eq. (85) 
was obtained from functional derivative methods that 
lead to results that satisfy conservation laws, 59 usually 
this method is applied to perturbative one-particle self- 
consistent schemes. In the TPSC equations, all the func- 
tional dependence on vector potential is in G' 1 ' . One may 
question whether this preserves conservation laws. The 
full Ward identity is derived in Appendix B where one 
also finds comments on why it cannot be used to find the 
vertex corrections in the limiting case we are interested 
in. Given the difficulty of computing the current-current 
correlation function alone for only the wave vector q = 0, 
it will be clear why the full Ward identity cannot be ver- 
ified. As a test of particle conservation, we focus instead 
on how accurately the /-sum rule, Eq.(Al9) derived in 
Appendix A, 

Xj.j m (*0n =°) = JjJ2 ( nk *) ' ( 87 ) 

ko- x 

is satisfied numerically. In this equation, the occupa- 
tion probability (rtk<r) is computed with G^ 2 \ while the 
left-hand side is the zero-Matsubara frequency current- 
current correlation function Xj^jx? Eq.(85), obtained 
with the functional derivative approach that gives us 
dX^ 2 -* / dG^ as irreducible vertex. 

Typically, using 500i as the cutoff Matsubara fre- 
quency, i.e. about 60 times the bandwidth, the above 



B. Optical conductivity 

The optical conductivity and the effect of the vertex 
corrections are different depending on which side of the 
quantum critical point the system is and on the tem- 
perature. On the left-hand side of the critical point, 
the conductivity changes qualitatively as the tempera- 
ture approaches the crossover temperature to the renor- 
malized classical regime T* . Here T* is defined as the 
temperature at which £ sp — £ t /j. As shown in Fig. 3(a), 
at high temperature it has a Lorentzian-like shape at low 
frequency, whether we look at the bubble alone, the bub- 
ble with the first vertex correction or with both vertex 
corrections. In this region of the phase diagram, the low 
frequency conductivity is smaller if the vertex corrections 
are included. Then, as seen in the left part of Fig. 3(b), 
as T is lowered toward T* but still above it, the effect 
of the first vertex correction is to strongly decrease the 
low frequency conductivity, while the second correction 
does not just compensate this effect, but make the total 
even higher than the bubble alone. At higher frequen- 
cies, in the right part of Fig. 3(b), a hump, not present 
without vertex corrections, appears in the total conduc- 
tivity. Note that at this temperature, the spin fluctuation 
frequency u sp is about T/2. When T is around T* , in 
Fig. 3(c), the hump is more pronounced and finally, when 
T < T* , in Fig. 3(d), it becomes a very distinct peak. At 
this temperature, a hump similar to the one seen with 
vertex corrections at higher temperature appears in the 
bubble term alone. At this temperature, T w 600uj sp and 
£ sp = 147 ~ 6£t/j. If we were to compare our result to 
experiments, assuming an energy scale t = 350meV, the 
hump seen around 0.5t would correspond to the feature 
observed in the mid-infrared frequency range in the op- 
tical conductivity of electron-doped cuprates. 79 ' 80 Note 
that another clear hump appears in the conductivity 
with vertex corrections in Fig. 3(d). At this doping and 
temperature, the spin fluctuations are incommensurate, 
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FIG. 3. Optical conductivity with and without vertex cor- 
rections at different dopings and temperatures. From (a) 
to (d) the doping is p = 0.17 < p c and the tempera- 
tures are (a) T = 0.3t > T* , (b) T = 0.08t > T* , (c) 
T = 0.06 w T* and (d) T = 0.04* < T* . The other dop- 
ings are p — 0.205 = p c at temperatures (e) T = 0.2t and 
(f) T = O.Olt and p = 0.32 > p c at (g) T = 0.2t and (h) 
T = 0.2t. All panels on the right-hand side contain a blow up 
of the low-frequency region of the rightmost plot. The sym- 
bols are only shown for a small fraction of the total number 
of points in the grids. 



hence their effect on the spectral weight at finite energy 
can be quite complex. This additional structure in the 
conductivity may be a consequence of this incommensu- 
rability. 

At the critical doping, Figs. 3(e) and 3(f), the low 
frequency conductivity is lower than the bubble result 
when only the first vertex correction is taken into ac- 
count, while it is higher with both corrections. This ef- 
fect is much more pronounced at low temperature in Fig. 
3(f). New secondary peaks also appear at low tempera- 
ture. Those peaks are at frequencies considerably smaller 
than the peak seen at p = 0.17 below T* and, as will be 
clear from the DC resistivity results in the next subsec- 
tion, there is no pseudogap regime at this doping. Those 



peaks may nevertheless be caused by correlations that 
are present far beyond the critical doping, at least at fi- 
nite temperature. The effect of those correlations are also 
clear in the resistivity results in the following subsection. 

Finally, at a doping higher than the critical doping, in 
Figs. 3(g) and 3(h), the conductivity with the first ver- 
tex correction can be almost the same as the bubble re- 
sult at low frequency both at high and low temperature, 
although this correction does not vanish in Matsubara 

frequency (at T = O.Olt, xZ)« (°) is about 7% of Xj°X (°) 
and x^ c } x (i27rT) is about 14% of x^^ttT)). Adding 
the other correction makes the low frequency conductiv- 
ity increase substantially. However, while the finite fre- 
quency conductivity stays finite when the doping is closer 
to p c or smaller, at low temperature, in Fig. 3(h), it van- 
ishes completely just after the main low frequency peak 
when both vertex corrections are included. This peak is 
also sharper and higher at high doping (not shown). As 
will be confirmed in the next subsection, this means that 
the system becomes closer to a Fermi liquid, although it 
has not yet reached this regime at this doping. When 
it does, the DC conductivity will be inversely propor- 
tional to T 2 at low temperature so that, to conserve the 
weight which is roughly constant with respect to temper- 
ature, i.e. to respect the /-sum-rule, the width has to be 
proportional to T 2 . At zero temperature, since there is 
no impurity scattering in our model, the low frequency 
part becomes a delta function. Note that an absoption 
band remains at finite frequency at the highest doping 
since there is still an incoherent part in the one particle 
spectrum. However, one must also note that the highest 
values of the conductivity in that band are about a thou- 
sand times smaller than the DC value and that the total 
optical weight of this band is about ten times smaller 
than the weight of the low frequency peak. 

The frequency region where the effect of vertex cor- 
rections is important becomes smaller as the doping in- 
creases. For p = 0.17, the difference between the different 
results vanishes around uj — At, at the critical point, this 
happens around uj = It and at p = 0.32, around u) = 1.5t. 

As the doping increases, the conductivity at low tem- 
perature becomes extremely sharp and it becomes very 
hard to do the maximum entropy analytical continua- 
tion. That is because a very fine grid must be used at 
low frequency, while one still needs a cutoff larger that 
the bandwidth. This makes the number of points in the 
real frequency grid explode, as well as the time for the 
optimization process. 



C. DC resistivity as a function of temperature and 
doping close to the quantum critical point 

Figure 4 shows five interesting doping regimes for the 
DC resistivity. The left vertical axis is in units of the 
Ioffe-Regel maximum metallic resistivity. The right ver- 
tical axis translates the result in micro-ohm-centimeters 
by taking d = 5 A as the interplane lattice constant. For 
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FIG. 4. Resistivity as a function of temperature for a) p = 
0.15 b) p = 0.17 c) p = 0.205 (critical doping) d) p = 0.26 and 
c) p — 0.32. The dashed vertical line in (a) and (b) indicates 
the temperature at which the internal accuracy starts failling 
as T decreases. Dashed curves in (c), (d) and (e) are the 
results of fits of the low T resistivity to the form AT + BT 2 . 



the temperature scale, we use t = 350meV. The resis- 
tivity without vertex corrections (bubble) appears as red 
circles, with the first correction, namely the second term 
in expression (85) or the Maki- Thompson-like diagram in 
figure 1, in purple squares (bubble-)- VC1) and the total 
resistivity appears as blue triangles. The insets of the 
last three figures display the low temperature behavior. 
The position of each doping is indicated by an arrow for 
each figure on a schematic "phase diagrams" with the 
quantum critical point and the crossover line. 

At the smallest doping, p — 0.15 in Fig. 4(a), one ob- 
serves at high temperature the expected Ioffe-Regel max- 
imum metallic resistivity saturation. The value is about 
2hd/e 2 . The resistivity is higher when vertex corrections 
are included at high temperature. Without vertex cor- 
rections, it increases with decreasing temperature below 
the crossover temperature to the renormalized classical 
regime T* , which is about 4Q0K (O.lt) at this doping. 
The effect of vertex corrections is dramatic at low tem- 
perature, essentially changing the resistivity from insu- 
lating to metallic. An important point is that this effect 
can only be obtain when both corrections are included. 
When only the first correction is included, the scattering 
effect of magnetic correlations is largely overestimated 
and the resistive behavior thus amplified. At p = 0.17, 
closer to the quantum critical point, Fig. 4(b) exhibits 
essentially the same behavior except that, given the over- 
all smaller resistivity, the saturation at high temperature 
occurs beyond the range displayed. At low temperature 
in Figs. 4(a) and (b), the resistivity with vertex correc- 
tions seems to extrapolate to negative values. However, 
as indicated by the vertical dashed lines in those figures, 
this happens in a region where the TPSC result is no 
longer quantitatively reliable because the internal accu- 
racy check starts to fail when the magnetic correlation 
length becomes too large compared to the thermal De 
Broglie wavelength. The results at those temperatures 
and dopings, namely, deep in the peudogap regime, can 
thus at most be considered as qualitative tendencies. 

At the quantum critical point, p = 0.205, one observes 
in Fig. 4(c) that, without vertex corrections, the resistiv- 
ity is quite linear at low temperature, as found previously 
in spin fluctuation theories. 21,81 When all the corrections 
are included, the most obvious effect is that the resis- 
tivity decreases at all temperatures. The linear behavior 
remains, at low temperature, but the vertex corrections 
tend to reduce the linear contribution and a T 2 behavior 
appears at a lower temperature than without vertex cor- 
rections. Note that, with only the first correction, there is 
a change in curvature. Also, the resistivity is larger than 
the bubble result at all temperatures instead of smaller 
as found when all vertex corrections are included. 

As the doping becomes larger than the quantum crit- 
ical doping, figures 4(d) and 4(e) show that a linear 
T behavior is still present at the lowest temperatures, 
but gradually disappears as the doping increases and the 
Fermi liquid- like T 2 behavior becomes dominant. The 
resistivity with the first vertex correction only is omitted 
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FIG. 5. Coefficients a) A and b) B in the fit of the form AT + 
BT 2 to the resistivity, as a function of doping, starting at the 
QCP. The inset in (a) shows the superconducting transition 
temperature as a function of doping estimated with the TPSC 
approach, from Ref. 82. The doping region relevant for our 
fits is indicated with dashed lines in this inset. 

in figure 4(e) because it is almost equal to the bubble 
result at all temperatures. 

The result of fits of the temperature dependence of 
the resistivity to the functional form AT + BT 2 over the 
range 0.008t < T < 0.05t (30K < T < 200K) is illus- 
trated in Figs. 5a and 5b. Those fits were done using 
two other resistivity curves in addition to those shown 
in figure 4, namely, at p = 0.22 and p — 0.24. The 
linear coefficient A decreases as one moves away from 
the quantum critical point. This decrease is correlated 
with the superconducting transition temperature, shown 
in the inset, estimated from TPSC calculations of the 
d x 2_ y 2-wave susceptibility 82 . As for the coefficient B. it 
seems to have a rapid increase as we move away from 
the quantum critical point (QCP) and then to remain 
roughly constant as the doping increases. Note however 
that this coefficient is hard to obtain precisely close to 
the QCP. For example, if we change the number of points 
used in the fit, the value of B in this doping region can 
change by about 20%. This is because at low tempera- 
ture, the region of interest here, the quadratic term has 
a very small contribution compared to the linear term 
close to the QCP. Also, there is always some noise in the 
resistivity, a consequence of the fact that those results 
are obtained by analytical continuation of Matsubara re- 
sponse functions, a procedure very sensitive to finite pre- 
cision noise in these response functions. Therefore, the 



value of B for the two smallest dopings in figure 5(b) 
are rough estimations. However, as the doping increases, 
the quadratic contribution becomes more important and 
thus B is more accurate. 



IV. DISCUSSION 

Starting with the self-energy at the second level of ap- 
proximation, (72) the only non-straightforward part of 
the derivation is the neglect of the functional derivatives 
of the irreducibles vertices T sp and Y c h in the deriva- 
tive of E^ 2 ) with respect to the field, or with respect to 
G (1) . However, as was explained in subsection II D, fol- 
lowing equation (73), the vanishing of those contributions 
to the conductivity is exact. Therefore, starting from the 
approximation (26) for the self-energy that contains the 
local spin and charge irreducible vertices r sp and T ch , 
Eqs. (51) and (52), the rest of the calculation is exact. 

In fact, based purely on our numerical results, with 
the Green's function at the second level of approxima- 
tion and the irreducible vertex generated from the cor- 
responding self-energy, we would have strong reasons to 
suggest that the /-sum rule is satisfied exactly in our ap- 
proach. At all dopings and temperatures, that sum rule 
is satisfied with very high accuracy and the precision in- 
creases with the cutoff Matsubara frequency. Since the 
/-sum rule is a consequence of particle conservation (see 
appendix A), this means that our approach is consis- 
tent with particle number conservation. This result is a 
consequence of our use of functional derivatives to cal- 
culate the Green's function and correlation functions, as 
in conserving approximations. 59 However, unlike in those 
approaches, there is no self-consistency at the one parti- 
cle level, i.e., our self-energy £( 2 ) is calculated with 
and not G^. Thus one-particle self-consistency is not 
necessary to ensure particle conservation. 

The saturation of the resistivity to the Ioffe-Regel max- 
imum value at high temperature is a general constraint, 
this time based on general physical considerations, that 
is satisfied by our approach and not by Boltzmann-based 
approaches. At high doping the Fermi-liquid T 2 resistiv- 
ity is also recovered, as expected. 

From the fact that the first vertex correction, the Maki- 
Thompson-like diagram in Fig.l, is sufficient to respect 
the sum rule, one would be tempted to assume that the 
other type of vertex corrections will not contribute, as 
was done in previous calculations. 27 It is clear from our 
results that this is not the case. In the pseudogap regime, 
the Aslamasov-Larkin-like contribution has a drastic ef- 
fect on the DC resistivity, making the system metallic 
instead of insulating. If only the first correction is in- 
cluded, we would wrongly assume that the system is 
even more insulating than without any correction. On 
the right-hand side of the quantum critical point, it is 
also very important to include both vertex corrections 
since their total effect is to reduce the resistivity, while 
it would increase, for a large range of doping, if only the 
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first correction was taken into account. Also, when both 
corrections are included, the T 2 Fermi liquid resistivity 
is recovered at a lower doping. The importance of in- 
cluding both corrections is also very clear in the optical 
conductivity, especially on the left of the quantum criti- 
cal point, where the result is qualitatively very different 
with only the first vertex correction. 

One may ask whether it is enough to include only those 
two vertex corrections. To answer that question we recall 
that our expression for the current- current response func- 
tion (85) is not derived perturbatively, but using func- 
tional derivatives. Since the calculation of Xj x jx ^ s exact 
starting from the TPSC self-energy, Eq.(26), as long as 
the TPSC approach is valid, all the terms needed are the 
ones that we have used. The region of the phase diagram 
where the TPSC result breaks down is deep in the pseu- 
dogap regime, when some parts of the Fermi surface are 
destroyed and thus the one-particle spectrum is dramat- 
ically different from the non-interacting one. Therefore, 
in this regime, the TPSC results should be regarded as 
more qualitative than quantitative. 

Since we have not used realistic band parameters, we 
cannot directly compare our results with any real mate- 
rial. However, some tendencies can help us understand 
experiments and provide some hints on what would be 
interesting to investigate using more realistic band struc- 
ture. 

First, there is a pseudogap in the TPSC approach. For 
electron-doped cuprates, it was shown that the pseudo- 
gap has properties predicted by this approach, namely 
it appears when the antiferromagnetic correlation length 
becomes larger than the single particle De Broglie wave- 
length, or the mean free path. 55 The hump around 
co = 0.5t, seen in figure 3(b), (c) and (d), thus corre- 
sponds to the analogous feature seen, for example, in 
Nd 2 - x Ce x Cu04 in the mid-infrared energy range. 79,80 If 
we use t = 0.35eV, the location in energy of this structure 
is also in the correct energy range. Those optical con- 
ductivity results and previous comparison of TPSC spec- 
tral weights with photoemission experiments on electron- 
doped cuprates add to the evidence that those materials 
are well described by the accurate solutions of the Hub- 
bard model provided by the TPSC approach at interme- 
diate coupling. 44,83 

An interesting aspect of our results is the linear low- 
temperature resistivity observed at the quantum critical 
doping. As mentioned in section IIIC, linear T resis- 
tivity has already been obtained in spin-fluctuation the- 
ory within the so-called self-consistent renormalization 
approach. 21,81 This is discussed in the review Rcf.24. 
However, this theory is based on the variational approach 
to the Boltzmann equation. It has been shown with a 
better variational ansatz that the resistivity should be in 
T 2 . 22 That variational ansatz is better because it takes 
into account that hot regions on the Fermi surface that 
are strongly influenced by spin scattering should be short- 
circuited by the cold regions. Indeed, it is the conductiv- 
ities of the different Fermi surface regions that are inte- 



grated and not the resistivities. Thus, like parallel resis- 
tors, the total resistivity depends only on the more con- 
ductive regions. In our results, vertex corrections tend 
to make the resistivity more quadratic, but a linear term 
remains at low temperature. This result tells us that, at 
least for the parameters used in this work, the region of 
the Fermi surface with the lowest resistivity has linear 
T resitivity. This is coherent with our observation, from 
spectral density calculations (not shown), that the whole 
Fermi surface seems incoherent, or hot, over a range of 
dopings that extends beyond the QCP at the tempera- 
tures considered. 84 

Another important question is whether, in the quan- 
tum critical region above the QCP, the whole Fermi 
surface is always incoherent, as suggested by some 
experiments. 85,86 This is possible because what charac- 
terizes this region is that the magnetic fluctuations are 
strong, but not yet strong enough to destroy any part of 
the Fermi surface. In other words, instead of having spin 
excitations that are well defined, they are broad in energy 
and wave vector, which means that they affect large parts 
of the Fermi surface and possibly all of it. To understand 
the transport, it is important to make the distinction be- 
tween this regime and the pseudogap regime, where the 
wc^ll defined hot spots appear. This distinction is not 
made in Rcf.22. In the present work, the incoherence 
of the whole Fermi surface may come however from the 
fact that we have included only nearest-neighbor hop- 
pings and that large parts of the Fermi surface are al- 
most nested. More calculations of the conductivity using 
second and third nearest-neighbor hoppings will be nec- 
essary to verify if this global incoherence of the Fermi 
surface in the quantum critical region is universal. Note 
that this is a subtle question that cannot be answered 
using approaches with adjustable parameters, like those 
of Refs. 21 and 22, which can give qualitatively different 
results depending on how those parameters are chosen. 

Another interesting point about the linear term in 
the resistivity is its correlation with the superconduct- 
ing transition temperature T c . This correlation seems to 
be present in all the unconventional superconductors. It 
has been observed in the cuprates, the pnictides, and the 
organics 87-90 . In those materials, it seems that the linear 
coefficient A disappears exactly at the end of the super- 
conducting dome on the overdoped side. Which strongly 
suggests a common origin for linear resistivity and su- 
perconductivity. In the Hubbard model, in the weak to 
intermediate coupling regime, the TPSC method finds 
that superconductivity and linear resistivity also have 
the same origin, namely the interaction of quasiparti- 
cles with antiferromagnetic fluctuations. However, from 
Fig. 5, it does not seem that A disappears completely with 
T c , which vanishes at around p = 0.25 as shown in the 
inset of the figure. Note that there is some uncertainty 
in the values of A that could come from a small system- 
atic error in the resistivity. This is possible because it 
is difficult to obtain very precise analytically continued 
results. Nevertheless, it is clear that the most important 
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drop in A is before T c vanishes. 

One striking result we obtain in the pseudogap regime 
on the left side of the QCP is the change of the resistivity 
from an insulating to a metallic behavior when the vertex 
corrections are taken into account. This is an effect of the 
Aslamasov-Larkin-like contribution which, as mentioned 
at the end of subsection II D, gives a positive contribu- 
tion to the conductivity when the density of states below 
the Fermi level is larger than above. In the present case 
the Fermi level is just above the Van Hove singularity. 
The density of states is therefore extremely asymmetric 
at the Fermi level and that explains why the AL term is 
so large that it counters both the effect of the one-particle 
self-energy that makes the bubble result insulating and 
the effect of the Maki-Thompson-like term that tends to 
make the system even more insulating. From preliminary 
results with second and third nearest neighbor hoppings 
if and t" , we notice that the second vertex corrections 
cannot always compensate the first correction and the re- 
sistivity for p < p c can be higher than the bubble result 
even with both vertex corrections. This happens when 
the Fermi level is farther from the Van Hove singular- 
ity, so that the density of states is much less asymmetric 
around w = and thus the AL term is weaker. Those 
results tell us that the behavior observed in the pseudo- 
gap regime on the left side of the QCP in Figs. 4(a) and 
4(b) is not universal. However, on the right-hand side 
of the QCP, the results should be qualitatively the same 
since eventually Fermi-liquid physics dominate. That is 
effectively what we observe so far. 

Another interesting question is the effect of disorder 
on the importance of vertex corrections. In our case, as 
mentioned in section II D, if the vertices do not depend 
on wave vector, i.e. they are isotropic, they have no ef- 
fect on the conductivity. Adding disorder would probably 
make the vertices more isotropic, and therefore their ef- 
fect should be smaller. It would be very interesting to see 
if, by adding enough disorder, we could reduce the effect 
of vertex corrections to the point where the resistivity 
recovers its insulating behavior in the pseudogap region. 
This would also provide one possible explanation for the 
fact that, in less clean cuprates, such as L^-zSr^C^C^ 
the resistivity increases below T* , while it increases in 
cleaner systems such as YBa2Cu307_ y . That contrast- 
ing behavior can be seen in Ref. 91 when, for the latter 
compound, we take the pseudogap line to end at opti- 
mal doping. However, from our preliminary results with 
hoppings t' and t" relevant for cuprates, we have noticed 
that, on the hole-doped side, the qualitative behavior 
with and without vertex corrections can be inverted with 
respect to the results of Figs. 4(a) and 4(b). That is, 
the resistivity without vertex corrections decreases be- 
low T* , while it increases with vertex corrections. By 
changing the band parameters and the doping we can 
therefore change dramatically the transport properties 
in the pseudogap regime. 



V. CONCLUSION 

To satisfy current conservation, conductivity calcula- 
tions must include vertex corrections that are consis- 
tent with the self-energy. A systematic way of achiev- 
ing this proceeds with functional derivatives with respect 
to the vector potential. We have shown how this ap- 
proach can be generalized to the non-perturbative TPSC 
approach. The various terms of the resulting algebraic 
expression, Eq.(85), have the physical interpretation il- 
lustrated schematically in Fig.l. One type of vertex cor- 
rection has the structure of the Maki-Thompson term 
in fluctuation superconductivity, while the other has the 
structure of the Aslamasov-Larkin term. These diagrams 
contain many elements that are not computed pertur- 
batively, as for example the irreducible spin and charge 
vertices, and the vertex corrections. With this approach, 
the /-sum-rule is satisfied in principle exactly. We veri- 
fied the agreement to the accuracy of the numerical cal- 
culations (a typical relative precision of 1CT 7 ). 

The numerical evaluation of the conductivity with ver- 
tex corrections can be done only if FFT and other ad- 
vanced numerical algorithms are employed. Brute force 
calculations arc impossible with any kind of computing 
resource presently available. Analytical continuation is 
performed with a specialized maximum entropy tech- 
nique that we have described in detail in appendices, 
along with all other algorithms. 

We have shown that our approach allows us to compute 
the optical conductivity and DC resistivity of the nearest- 
neighbor two-dimensional one-band Hubbard model in 
a variety of regimes, without adjustable parameter or 
phenomenological assumption. There is no need to as- 
sume the existence of quasiparticles, as is the case in 
Boltzmann equation approaches. 

For illustrative purposes, we have presented the results 
of calculations for U = 6t and nearest neighbor hoppings 
only. For the DC resistivity, we find at high tempera- 
ture that it saturates at the Ioffe-Regel value, namely 
when the mean-free path is of the order of the inter lat- 
tice spacing. The existence of this Ioffe-Regel limit is 
usually assumed on phenomenological grounds. Here we 
have demonstrated it. That limit may be exceeded in 
strong coupling 92,93 , but that is beyond the TPSC ap- 
proach. Vertex corrections have a dramatic influence for 
dopings smaller than the quantum critical doping. They 
can change the temperature dependence of the resistiv- 
ity from insulating to metallic when the antifcrromag- 
netic correlation length becomes larger than the thermal 
de Broglie wave length, namely in the pseudogap regime. 
At the quantum critical point, the resistivity is linear 
at low temperature, although vertex corrections tend to 
reduce the temperature range of the linearity compared 
to the calculation with the bubble only. At low tem- 
perature, the linear term persists at dopings larger than 
the quantum critical point, although we cannot exclude 
that it disappears at temperatures lower than those ac- 
cessible to us. The coefficient of this linear term is also 
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correlated with the vanishing of the superconducting T c 
obtained with the TPSC method, in qualitative agree- 
ment with experimental results on the cuprates, the pnic- 
tides and the organics. 88,89 For dopings equal or greater 
than the critical doping, the resistivity with all vertex 
corrections is always smaller than the simple bubble re- 
sult. In general, for most of the dopings and tempera- 
tures considered, the first vertex correction has the ten- 
dency to increase the resistivity, while the second has 
the opposite effect. Therefore the results of each type 
of term are always very different and one cannot neglect 
the Aslamasov-Larkin-like contribution in those regimes. 
The latter contribution vanishes only when there is exact 
particle-hole symmetry. 

We observe in the optical conductivity that vertex cor- 
rections are important at all the dopings considered and 
that no term can be neglected. The effect is the strongest 
for dopings smaller than the critical doping near and be- 
low the pseudogap temperature. We also observe that the 
frequency at which the results with and without vertex 
corrections cease to differ decrease with increasing dop- 
ing. The hump structure in the mid-infrared frequency 
range, related to the pseudogap, is observed both with 
and without vertex corrections, but with different ampli- 
tudes at a given temperature. At the quantum critical 
point and beyond, the effect is important both at high 
and low temperature though more important at low tem- 
perature. The low frequency part of the conductivity is 
higher with the vertex corrections and increases quite 
rapidly with doping. 

This work is presently being extended in several direc- 
tions. For example, one can study more realistic band pa- 
rameters for high-temperature superconductors and in- 
vestigate the connection between the single-particle scat- 
tering rate along the Fermi surface and the temperature 
dependence of the resistivity. The sensitivity of the re- 
sistivity to the details of the model in the pseudogap 
regime would also be interesting to investigate. Using a 
similar approach, one can also envisage calculating the 
thermopower and other transport properties. 
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Appendix A: /-sum rule for the conductivity 

In this appendix, we derive the general expression for 
the real part of the conductivity, Eq. (19), and we show 
that the zero-Matsubara frequency value of the current- 
current correlation function suffices to check numerically 
the validity of the /-sum rule. We begin from the conti- 
nuity equation, 



dp(r,t) 
dt 



+ V-j(r,i) =0, 



which, in Fourier space, reads 

-wp(q,w) + q j(q,w) =0. 
If j(q,w) =j x (<l,u)x, we have 

Qx3x(<l,u) = wp(q,w) . 



(Al) 

(A2) 
(A3) 



Using space translational invariance, the two-particle 
spectral function corresponding to an observable A can 
be written formally as 

xl 1 (q,w) = -^ f <[A(q, W ),A(-q,- W )]). (A4) 

with T — > oo, so that, from the continuity equation 
Eq. (A3) we have the relation between current and charge 
correlation functions 



(A5) 



and thus, 



J v— — = gj v wx - (q ' w) - (A6) 

The right-hand side can be obtained from equal-time 
commutators since 



dw / . d 

— w Xpp(q^) = I *■ 



dw 



dt, 2ir 



e-*" l 2 X " (q,w) 



and, by definition, 



t=o 
(A7) 



/ 



duj 



e-^x'^uj) = -([p(q,t), P (-q,0)]), (A8) 



2tt 
so that 



/ 



w 



11 d 

1 1 



t=0 



([[p(q),if](i),p(-q,0)]> 



t=o 
(A9) 



Taking for H the Hubbard hamiltonian and calculating 
the commutators, we get 



/ 



dw x'jjM^) = 11 
" <ll N 



LO 



^(ek+q + £k- q - 2e k ) (n kcr ) . 

(MO) 



krx 
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We are interested in the long wave length limit, namely 
q — > 0, so that 

, <9e k <9e k 

€k+ ^ ek + W x qx + W y qy 

l<9 2 e k 2 19 2 e k 2 d 2 e k 

'/.'/.,/ (All) 



+ n-579-E+o-S79-K + 



2 dki Hx 2 afc 2 dk x dk v 



and, if we consider the longitudinal conductivity q = q x x, 
we obtain 



f doj \ ; 1 \- 9 2 e k , . 

J T « = jvf a^ (nkCT> 



kcr 



(A12) 



using the earlier definition Eq. (12). 

To derive the general results for the real part of the 
conductivity Eq. (19), we begin with the spectral form 
for Xj x j x (lx,u) that reads 



doj' x'LiAQx'V) 



■ oj — irj 

so that, from Eq.(18) and Eq.(Al2), we have 

( k x) + Xj x ja> M 



(A13) 



&xx(qx, oj) 



i(w + irj) 



i(oj + irj) 



TT OJ' 



+ 



r duj' x'j 3 Mx^') \ 

J TT oj' — oj — IT] ) 



■ OJ — ITj 

1 f doj' (u + iri)x'- xi (q x ,uj') 



i{uj + irj) J tt oj' (oj' — oj — irj) 

1 fd^ xjjJ^O 
i J 7r oj'(oj' — oj — irj) 



Substituting in the form found earlier Eq.(A12), the /- 
sum rule for the conductivity is 

fit). , , 1 9 2 e k . . 
7 _ -Reo xx { qx ,oo) = -g w (n kff ) , (A17) 

for g x — > 0. Since the spectral form for Xj x jAlx,i<ln) is 

^.fe,*)=r^^%^. (Am 

J-oo T oj -iq n 
we have the desired result, 

Xj.j. (fe, «9n - 0) = — -gjjr ( n k ff ) • (A19) 

kcr x 

It is thus sufficient to look at the zero Matsubara fre- 
quency value of Xjxjx to check whether the sum rule is 
satisfied (assuming that (k x ) has been calculated). 

Since the results Eq.(A17), or equivalently Eq.(A19), 
are a consequence of the continuity equation, it means 
this sum rule must be respected when there is conserva- 
tion of particle number. 



Appendix B: Ward Idendity 

In this appendix, we derive the general Ward iden- 
tity that follows from charge conservation. It is far 
too complicated to be verified in full generality numer- 
ically within our approach. We also indicate that the 
Ward identity suffices to find the vertex correction sim- 
ply in cases where frequency and momentum variations 
are smooth. These assumptions are not fulfilled in our 
case hence these simplifications cannot be used. 

We begin from linear response theory in imaginary 
time. We omit the diamagnetic term, which is not rele- 
vant for the present discussion. We find, 

(A14) (jq (r)) = [ dr' (T T j q (r) j_ q (/)) • A q (/) (Bl) 



and since 



' — +mS{J-u), (A15) 



oj' — uj — irj oj' — oj 



we obtain the desired result 



D / \ X'j x j x (<lx,u) 

Rea xx {q x ,oj) = 



■oj 



(A16) 



where j is the paramagnetic current in the interaction 
representation. Take a single (bosonic) Matsubara fre- 
quency for the vector potential 



A q (r')=TA q (g m ) e -^ 



(B2) 



and extract the corresponding Matsubara frequency for 
the current, then the quantity to evaluate is 
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Gq(9m)> =T J* dr j* dr'e^-r') (T T j q (r)j_ q (/)) • A q (q m ) 

=T f dr f dr'e^-r') 1 £ £ Vk£k (B3) 

C Lr ( T ) C k+qCT (t) C k'+qcr ( T ') C k'cr ( T ')^ ^k'^k' ' A q (<?„,) 

where, in the last equality, we have also used spin conservation with the fact that only the connected piece will 
contribute because the average current in equilibrium vanishes. 

The Ward identity that we need can be obtained from the single spin component version of following equality that 
can be derived from current conservation, 94 

^r + ( e k+ q -£k) (TrcLWck+ q ,( r )4'+ q a( T i)ckv(r 2 )) ^ 

= S (T - Tl) Gk'cr (r 2 ~t)-5{t - T 2 ) Gk'+qo- (r - n) . 

Current conservation, that we saw in the first two equations of appendix A, would give a vanishing right-hand side 
were it not for the theta functions whose derivatives give delta functions and ultimately equal-time commutators that 
can be evaluated. Define the fermionic Matsubara frequency four-point correlation function 

K (k,k + q;fc m ,k' + q;*4,k') 

EE £ dn J* d T 2 e ik ^-Tl) e ik' n (T 2 -r) (^J^ (t) Ck+ ^ (r) c t ;+qCT (n) Ck , ff . 

Then, taking the same fermionic components of the Ward identity, integrating by parts, it takes the form 

IM fc ™ + ik n) + ( £ k+ q - £k)] K (k, k + q;fc m , k' + q;k' n , k') = G k - ff (k' n ) - G k < +q(J (k m ) . ( B6 ) 

k 

We need the amputated function that is summed over all wave vectors to compute the current- current correlation 
function. So let us define the charge and current three point vertices (for a single spin component), valid in the small 
q limit 

T p (k m , k' + d;k' n , k') = - ^ A CT (k, k + q;fc m , k' + q;^, k') Gj a (k' n ) G k / +qff (k m ) (B7) 

k 

q • Tj (k m ,k' + q;k' n , k') = -q- £ V k e k A ff (k, k + q;fc m , k' + q;k' n , k') G k / CT G k / +qff (fc m ) . (B8) 

k 

Then, in the long wave length limit, the general Ward identity Eq.(B6) can be rewritten for the three-point functions 
as follows 

(zfc m - ik'J r p - q ■ Tj =G k / +q(T (fc m ) - G k / ff (tfj 

= («fc m - ik' n ) - q • V k /£ k ' - S k ' +qCT (fc m ) + S k v (fcjj 

We can obtain the four-point function that we need in the expression for the current Eq.(B3) from the fermionic 
Matsubara expression Eq.(B5) for the four point function at n = r 2 = r' as follows 

(j q (<U) =t /V f^T'e i '»M'-VVrVrVe- il »Me- it »( / - r ) 

J ° 7 ° ^ k, (7 k' fc m fc^ 

x V k e k A ff (k, k + q;fc ro , k' + q;k' n , k') V k ,e k < • A q (g m ) (BIO) 
= ^EE T E Vk£kA - (k»k + q;fc' n + « m ,k' + q;tf n ,k')V k /e k / ■ A q (<z m ) . 

k,<7 k' k' 



The sum over wave vectors k allows one to rewrite the latter in terms of the current vertex defined in Eq.(B8) 

2 

'N 



(j q =o (q m )) = ~ E T E r ^( fc « + «™ k '; fc - k ') Gk ' Gk ' + *») Vk ' £k ' • A i=° («™) • ( BU ) 

k' k' 



We have performed the sum over spins, which explains the factor two, and taken the q = limit first to represent a 
constant electric field. 

I 

There is a simple case where the Ward identity suf- ficcs to find the vertex correction. Consider the Ward 
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identity for the three point function Eq.(B9). The 
rest of this Appendix is correct only if we can assume 
that in the q — )■ finite frequency ik m — ik' n limit, the 
charge vertex T p (k m ,k ; + q;k' n ,k') and current vertex 
Tj (k m ,k' + q;k' n7 k') are analytical in q and have a fi- 
nite limit at q = 0. It can be checked that this is the case 
for the non-interacting system. In that case, all the non- 
analyticities are contained in the product of Green's func- 
tions appearing in our last expression for (j q= o (<7m)) ■ 
While these kinds of analytic properties can be assumed 
for Fermi liquids, this is not appropriate in our case, that 
is more singular in the presence of strong antifcrromag- 
netic fluctuations. 

In the case where we can assume analyticity, the expan- 
sion for T p in powers of q must begin at order q 2 because 
it is a scalar. Then, the long wave length limit of the ver- 
tices can be found by identifying the coefficients of the 
(ik m — ik' n ) and of the q on the left and right-hand side 
of the Ward identity for the three point function Eq.(B9). 
For the charge three point function we thus find 

r p ( km , ko = i - Ew : grj {K) . (Bi2) 

while for the current vertex, taking the q — > limit, one 
obtains 



Tj (k m ,-k';k' n ,-k') = V k <e k , + V k >5W {k ri 



(B13) 



The first term in the last two equations is the bare vertex 
and the last term the vertex correction. These results for 
the vertex corrections are valid, for example, for impurity 
scattering or electron-phonon interactions where the as- 
sumptions upon which they were derived are valid, 95 but 
not in our case where the gradient of the self-energy with 
respect to wave vector can be very large, for example 
near hot spots. The above two equations by themselves 
are often called the Ward identities. 95 

Note that, at zero external Matsubara frequency, as 
long as the charge vertex does not diverge, Eq. (B13) 
becomes exact for k m = k' n . This form for the vertex 
generates ladder diagrams only, demonstrating that the 
Aslamasov-Larkin contribution vanishes, as found nu- 
merically in Sec. Ill A. 



Appendix C: Fast-Fourier transforms, cubic splines 
and asymptotic expansions 

The use of fast- Fourier transforms (FFTs) is absolutely 
essential for this calculation. We begin by describing 
their use for obtaining the Lindhard function, Eq.(22), 
and the self-energy, where the convolutions that render 
the use of FFTs possible are apparent. The convolutions 
are not so apparent in the case of the conductivity, par- 
ticularly the Aslamasov-Larkin like terms, that require 
an elaborate discussion in the third subsection. 



1. The Lindhard function 

The Lindhard function Eq.(22) is needed to compute 
the spin and charge susceptibilities Eq.(20) and Eq.(21), 
which is defined in a more explicit way as 



with 



(q,*g„) = -2- J2 G (1) (k+q,zfc„ + z 9n )G (1) (k, J fc„), 



k,ik„ 



G«(k,ifc„) = 



1 



(CI) 



(C2) 



where e k is the bare particle dispersion relation, k n = 
(2n + 1)ttT is a Fermionic Matsubara frequency and fio 
is the bare chemical potential. Expression (CI) could be 
calculated by first performing analytically the sum over 
Matsubara frequencies to obtain 



Xo(q, iq n ) 



2 /(&) ~ /(&+«,) 
N ^ iq n + Ck - Ck+q ' 



(C3) 



where £ k = e k — fi an d /(£k) is the Fcrmi-Dirac distri- 
bution, and then integrating numerically over k for each 
value of (q, iq n ). While it may seem that we have saved 
some work by doing exactly the Matsubara sum, by do- 
ing so we do not use a property of Eq.(Cl) that can make 
our calculation much easier. This is the fact that this ex- 
pression is a convolution, and convolutions can be calcu- 
lated in a very efficient way using fast-Fourier transforms 
(FFT). For instance, Eq.(Cl) can be written as 

xo(q, iq n ) = 

- 2 / dTe^ T ^e- jq ' r G (1) (r,T)G (1) (-r,-T), (C4) 
Jo 

where j3 = T~ x . 

To compute Xo(q, *9n) on a grid of size NqN qri us- 
ing the form where Matsubara sums have been done, 
Eq.(C3), one needs to do this number of integrals over 
the Brillouin zone, a number of operations that scales like 
N 2 Nq n if we consider that each integral scales like iV q . 
On the other hand, using the convolution form Eq.(C4), 
we need to do a two-dimensional FFT on G^^(k, r) 
to obtain G^\t,t) (here, G^^k, r) is known analyti- 
cally), a task that scales like TVqiV^ log-/V q , and then 
the three-dimensional FFT in Eq.(C4) that scales like 
N^N^ log NqN qn . We therefore have a gain proportional 

tO N J log NqNq„. 

There is one delicate point. Eq.(C4) contains a con- 
tinuous Fourier transform on the imaginary time r while 
FFTs are discrete transforms. The simplest thing to do 
would be to dicretize t and perform an ordinary FFT. 
This would give acceptable results for the low Matsubara 
frequencies but the high frequencies would be completely 
wrong since Xo(Q:*9n) would be periodic in q n , while it 
has to decrease as a series in even powers of l/q n . To 
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solve this problem we make use of cubic splines to ap- 
proximate the integrand between the discrete imaginary 
time points and perform a continuous Fourier transform 
on this spline. In appendix E, we show how to compute 
the continuous Fourier transform of a cubic spline using 
in fact only a discrete Fourier transform. Let us consider 
first the imaginary time Fourier transform in Eq.(C4). 
Using the form (E7) given in appendix E, we obtain 



Xo(r,T = P) -Xo(r, 



0) 



o 2 

'in 



N-l 



+■ - 5>i+i( r ) ei, " T< - (C5) 

3=0 



Qn 



where N is the number of intervals in the imaginary time 
grid, At is the size of an interval, Sj(r,r) is the cubic 
polynomial in the j th interval and (r) is the third 
derivative of Sj(r,r). 

Notice that Eq.(C5) contains the discrete Fourier 
transform of (r) , so that the result of this transform 
itself will be periodic in q n . But the factor l/(q n ) in 
front makes this term important only for low frequen- 
cies. Expression (C5) contains the derivatives of Xo( r )T") 
with respect to r at the boundaries. It is explained in 
appendix E how those derivatives are used to complete 
the linear system that must be solved to obtain the spline 
coefficients. They thus have to be calculated before the 
spline and must be known when Eq.(C5) is used. Here 
since xo( r , T ) = —2G^ (r, t)G^\— r, — t), and we have 
analytical expressions for G^\k, r), it is straightforward 
to calculate this derivative. For < r < /3, we have 



G< 1 >(k,T) = -e- £kT /(-&) 
= - e ^-)/(£ k ) , 

so that 

5G«(r,r) _ 1 ^ _ 



(C6) 



dr 



N ^ 

k 



(C7) 



>/(&)■ 



Then, still for < r < P, 

G«(k,-r) = e^/(a) 

and, therefore, 



(C8) 



9GW(-r,-r) 
dr 



^5' 



r ae« kT /(a) • (C9) 



If we have inversion symmetry, 

G (1 \-r,-r) = G^(r,-r). 

then 

dGW(-r,-T) 



(CIO) 



dr 



^E eikr ^ kT /(e k ), (en) 



and we have 

dGW(-T,-T) 



flGW(r,^-T) 

dr dr 



2. The self-energy 



(C12) 



The next function we can calculate using FFTs is the 
self-energy Eq.(26) that can be written as 

4 2) (Mfc„) = Un_ a 

+ J d T e tk ~ T e-^ r V{-r, -t)G^ (r, r) , (C13) 



where 



U 



V(v, t) = - [3U spXsp (v, t) + U chXch (r, r)] . (C14) 

To calculate V(r, t) accurately, we can use the fact 
that x S p(q, iq n ) and x c ft(q, iq n ) approach assymptoti- 
cally Xo(q, iqn) as q n increases and that Xo( r ; t) is known 
once Gi 1 ' (r, r) is. Then V (r, r) is computed from 



V(t,t) = 



N 



3U sp Xsp(q,iq n ) + U ch Xch(q,iqn 

- {3U sp + U ch ) xo(<i,iqn) 

+ ^(3U sp + U ch ) X o(r,T), (C15) 

where xo( r > T ) — —^^(rj^G^^—r^—T). Because the 
asymptotic part is removed in the Fourier transform, it 
converges as the transform of l/(q n ) instead of l/(g 2 ), 
so that a smaller cutoff can be used. 

In Eq.(C13), there is a continuous Fourier transform 
so that we have to use cubic splines to represent the r- 
dependence of the integrand. To compute the splines and 
then use formula (E7) we need to compute the derivatives 
of V(— r, — t)G^o\y, t) at r = and r = j3, i.e. those 
derivatives of V(— r, — r) and G^\t,t). The latter were 
already computed using Eq. (C7) to calculate xo- F° r 
V(— r, — t), since this function is symetric with respect 
to t = (3/2, we only need the derivative at t = 0. If we 
differentiate expression (C15) with respect to r and set 
t = 0, we notice that the sum disapears since Xsp(q ; *<Zn) 
and Xc/i( ( 1j iqn) are even functions of q n while the deriva- 
tive makes a factor q n appear in the sum. In other words, 
at r = (and thus r = /3), the derivatives of x SJJ (r), 
Xch( T ) and Xo( T ) are equal, which comes from the fact 
that they have the same asymptotic limit as can be seen 
from Eq.(E7) (with k = q n ), considering that there are 
no odd terms for those functions. Therefore, we have 



dV(r,T) 



dr 



U , sdxo(r,r) 
— (SUsp + U ch ) 



T = 



dr 



. (C16) 



T=0 
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Once the self-energy Eq.(C13) is obtained, we have to 
calculate the interacting chemical potential /i to define 
the Green's function G^ as 



G (2) (k,ifc„) 



1 



ik n - £k + A* - S ^ 2) ( k , *M 



(C17) 



For a given filling (n) , \i is defined implicitly by the equa- 
tion 



<rc>=2^G( 2 )(k,T = 0-) 



N 
T 



2 ^ £ e -ifc-o- G (2) (k) 



(C18) 



However, this expression is not very practical numerically 
because of the convergence factor e~ tkn0 that is not well 
defined numerically. We use instead 



2^ Yl [G (2) (kA)-G (1) (MM 



o. 



(C19) 



where it is assumed that G^ contains the chemical po- 
tential /iQ that gives the desired filling. No convergence 
factor is needed since the sum converges as the sum of 
V(fc 2 )- 



To use this formula we first have to compute G^ (k, r) 
from G^(k,ik n ) defined in Eq.(C17). To do that we 
cannot simply perform a Fourier transform with respect 
to ik n on G^(k,ik n ) including only a finite number of 
frequencies. This is because G^(k, r) has a discontinu- 
ity at t = that will produce oscillations (the Gibbs 
phenomenon) close to r = and t = /3 (remember that 
G(k, f3— t) = — G(k, — t)). This problem can be solved by 
using the asymptotic form of G^ 2 )(k, ik n ) to perform the 
transform on all frequencies. First, let us define G(k, r), 
the Fourier transform of some Green's function G(k, ik n ) 
on a finite set of frequencies, 



G(k,r,) = T^ ' e- ik "^G(k,ik n ) 



N T -1 



= e i*(N T -i)j/N TT e- l2 ^l N -G{\i,ik n _^). 

(C21) 



n=0 



3. The current- current correlation function 



Finally we come to the calculation of Xj^j^i Eq. (85). 
The first term, called the "bubble" contribution because 
of its bubble shape in figure 1, can be written as 

1 E f |r V E G(2) ( k > ^) G(2) ( k > ^ + 



N 







2 


I dre tg " T G (2) (k,r)G (2) (k,-r). 1 




\dk x J 








G^ f (k,ik n ) = 



where E means that the sum is finite, N T is the number 

of values of r and we have used k n = (2n + 1)ttT and 
tj = j/(N T T). The sum in second line has the form 
used in standard FFT routines. Now, to compute the 
self energy Eq. (C13), we have used Fourier transforms 
of cubic splines so that E( 2 )(k, ik n ) has the form Eq.(E7). 
The asymptotic form of G^ 2 ^(k, ik n ) is therefore 



(can) = 



e k - E in /(k,«fc„) 



gi I s 2(k) , s 3 (k) 
ik n ^ (jfe„) 2 " r (ifc„) 3 



(C22) 



where ek — £k — A*- The Fourier transform over ik n of this function can be done analytically using the residue theorem. 
For a function g(z) having only simple poles, we have 



^E« 



r g(ik m ) 



ik m 

which, applied to (C22), gives 

4 



J Z = Zj 

Ej Res [9(z)]f( Zj )e-^. 



< T < P, 

-B<t<0. 



C/(k,r) = Tle- z ' ,k) 7(T^(k))- 



fe(k) 



) - 2l (k)] . . . [z,(k) - ^(k)]^) - z, +1 (k)] . . . [ Zj (k) z 4 (k)] 
I 



(C23) 



(C24) 



where the minus sign is for < r < j3 and the plus polynomial 
sign, for —0 < r < 0. The 2j(k) are the roots of the 4 _ 



e k z 3 - siz 2 - s 2 (k)z - s 3 (k) . 



(C25) 
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Those roots are given by quite imposing, but analytical 
formulas. Finally, assuming that G^ ; f (k,T) is the finite 



Fourier transform of G- n ^-(k, ik n ) as defined in Eq.(C21), 
we use 

G( 2 '(k,r) = G( 2 )(k,r) + [G<£(k,r) - Gg(k,r) 

(C26) 

where the term between the brackets is the contribu- 
tion from frequencies beyond the cutoff used in Eq.(C21). 
This expression will be very accurate if the asymptotic 
behavior is well obeyed beyond the cutoff. Using the 
somewhat complicated expression (C22) as the asymp- 
totic form of G( 2 )(k, ik n ) may seem to complicate the 
calculations needlessly since the first term of the high 
frequency expansion of G^(k,ik n ) is identical to that 
of G^ 2 ' (k, ik n ). However, to compute the cubic spline of 
the integrand in Eq.(C20), and then its Fourier trans- 
form, we need the derivatives of G^(k, r) at r = and 
t = f3 and, for that purpose, we have to use a more ac- 
curate asymptotic form. The reason will be explained 
shortly. 

As we have just mentioned, we need to compute deriva- 
tives of Eq.(C26) with respect to t at the boundaries. Let 
us first rewrite expression (C26) using Eq.(C21), 

G( 2 >(k,r) =T^'e-^ |g< 2 > (k, ik n ) - Gg(k, ik n ) 



+ Gg(k,r) (C27) 



so that 



9G( 2 )(k,r) i 
dr I t=o 

r^(-i*n) [G (2) (k, i fc„) - Gg(k,ifc n ; 



+ 



9G^(k,r) 



dr 



(C28) 



where the last term is obtained from the derivative of 
Eq.(C24). If G^J f was taken to be the sum in (C27) 

would converge like the sum of l/(ik n ) 2 , but the sum in 
(C28) would not converge. Thus, G\ n f must have at least 
the same first two terms in its high frequency expansion 
as G^ 2 \ If we use Eq.(C22), the first five terms in the 
expansion are equal, the sum in (C27) thus converges as 
the sum of l/(ifc„) 6 , while the sum in (C28), as the sum 
of l/(ik n ) 5 . This gives us a very precise evaluation of 
G^ 2 )(k, t) and its derivatives. To obtain the derivatives 
at t = (3 we use the relation 



9G CT (k,T) 



<9r 



e k + Un- a - 



T=0 



8t 



(C29) 



r=0 



that is derived from the spectral representation of 
G(k,r). 



Once G^ 2 ' (k, r) and its derivatives are obtained, the 
integrals in expression (C20) are evaluated by comput- 
ing the cubic splines for G' 2 )(k, T)G^ 2 ^(k, — r) and us- 
ing formula (E7) for the Fourier transform. The rest of 
the calculation is simply a sum over the Brillouin zone, 
where, of course, it is preferable to use the symmetries of 
the system to save computational resources. 

The second term of Eq.(85) is 

V / k ± k 2 

xG ^(k 2 )G^(k 2 + l q n )^Jk 1 )^(k 2 ) 

x [iU sp Xs P {k 2 - fci) + U ch Xch(k 2 - fci)] 

(C30) 

where k\ + iq n = (ki, ik m + iq n ). If we define 
fn(k 2 ) = ^{k 2 )G {1 \k 2 )G^(k 2 + iq n ), 

V(k 2 - fci) = — [3U sp Xsp(k 2 - fci) + U ch Xch{k 2 - fci)] , 

(C31) 

then, Eq.(C30) becomes 

fel x 

x^/„(W2-fci). (C32) 

k 2 

The sum over k 2 being a convolution, we have 

Xl j Mn) = -lY,W x G{2) ^ Gi2) ^ + i ^ 
k 

xY^e- ik - l f n (l)V(l). (C33) 
i 

Note that the sum over 1 written in explicit form is 

f dre-^e^Ur^rMr^T). (C34) 

where the function /„(r,r) is given by 



xT^2 e- lk ™ T G^ (k, ik m )G^ (k, ik m + iq n ) . 

ik m 

(C35) 

bmce G«(k , ik m ) is the non-interacting Green's func- 
tion, the Fourier transform over ik m can be done analyt- 
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ically. Using (C23), we obtain, for q n ^ 0, 



T ]T e~ ikmT G^ (k, ik m )G^ (k, ik m + iq n ) 



ik n 



e iqnT _ I 



e"* T [0(r)(l-/(e k ))-0(-T)/(e k )] 



(C36) 



For g n = 0, there is a double pole so that the calculation 
with the residue theorem is slightly different. However, 
we can use the simple following trick, 



T £ e -*»* G« (k, ifem )» = T ^ e-*- ( Y 

ik \ 771 / 



~- (-e- e "* T [(l - /(e fc ))<?(r) - /(e~*)0(-r)]) 



e- ?kT r[(l - /(e k ))0(r) - /(e k )6»(-r)] + 



g/(gO 

<9e k 



(C37) 



For t > and q„ 7^ 0, we have 



i<?„ AT ^ dk x 



g n {r)h{v,T) . 



(C38) 



The functions /„(r, t) for different values of n are thus 
obtained by multiplying the n independent function 



(C39) 



by 3„(t) = (e 4<? " T — l)/(iq„). Using this result, we find 



= J rfr e '^9 n (T)^ e -' kl '^(r^)V(r jV 



(C40) 



Inserting this back into Eq.(C33), we obtain for the q n 7^ 
terms, 



T 
'.V 



4 E E w G(2) ( k ' ik ^ G{2) ( k ' ^» + 



= - — ^ E I? T E G(2) ( k > ifc ™) G(2) ( k > + 

«9n N ^ dk x ^ 

x ( / rf Te i ( ); ™+«»)^e- ik ' r '7i(r i ,T)y(r jl T) 

(C41) 

To make the convolutions more apparent, we define 
t7 (k, ik m ) — 



so that 



>(k,ifc ro ) I dre**»' T ^e- ik -^(r J - 1 r)V(r jl r) 



(C42) 



-if ^ E w ( T E G(2) ( k > J ( k > ^ + 

k 21 jfe m 

- G< 2 > (k, ik m + iq n )J(K ik m )) (C43) 
and we can use the convolution theorem to obtain 

XT xj Mn) = 

- — 4E ^ ( [ P dre-^G^(Kr)J(k,-r) 

iq n N ^ dk x vy 

- J dr e lq - T G {2) (k, r) J(k, -r)) . (C44) 

Note that, since all the Matsubara sums have been 
transformed into Fourier transforms that can be done 
with FFTs and that an FFT gives all N values of a N 
point transform at the same time, we obtain the values 
of Xj^j x (*9n) f° r an *9™ except iq n — at the same time 
when the last FFT is done. 

We still need to do the Fourier transform over r in 
Eq.(C42) using cubic splines. For that we need the 
derivatives of h(rj,T)V(rj,r) at r = and r = j3. The 
derivatives of V(rj,r) have already been used to com- 
pute the self-energy and are obtained using Eq.(C16). 
For h(Tj,r), the derivative is obtained from the defini- 
tion (C39). 
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Next, we have to compute J(k, — r). First, we define 
the function 

Q(k,ifc„) = / ' dTe ik « T Y i e- ik -''h(r j ,T)V(T j ,T) 

(C45) 

so that the definition (C42) reads J(k, ik m ) = 
(k,ik rn )Q(k,ik n ). The Fourier transform over r is 
done using the method described in appendix E. The 
asymptotic form of Q(k,ik n ) is 

n fir 4h \ , ga(k) , ga(k) , , 

Q inf (k, ik n ) = —— + —— + ——3 , (C46) 

so that, using Eq.(C22) for G?^(k, ik m ), the asymptotic 
form for J(k, ifc m ) is, 



1 



ik - f, - ( -ZL- 4- ■ S2(k) -I- 83 ( k ) ^ 
' gi (k) , 92 (k) , 93 (k) 



ifc n (z/c n ) 2 {ik n )^ 
and its Fourier transform is 

J m/ (k,-T)=£ e ^ k >VMk)) 

^ gi (k)[^(k)] 2 + 9 2 (k)^(k) + g3 (k) 



(C47) 



n^-[^(k)-^(k)] 



(C48) 



where the Zj(k) are the roots of the polynomial (C25). 
By analogy with our previous calculations, the function 
J(k, — r) is then obtained from 



J(k, -t) = e 4fc " T [J(k, zfc„) - J m/ (k, ifc„)] 

+ J m/ (k,-r) (C49) 



ik n 



where ^ means that the sum is done up to a cutoff 
frequency. 

Once these results are substituted in the expression for 
the Maki-Thompson term, Eq.(C44), we need the deriva- 
tives of G( 2 )(k,r)J(k, -r) at r = and r = /3. For 
G( 2 )(k,r), they are given by Eq.(C28) and Eq.(C29). As 
for the derivatives of J(k, — r), they are obtained by dif- 
ferentiating Eq. (C49) . 

Finally, we need to evaluate separately the q n = term 
for t > 0. From Eq.(C37), we have 



/o(r,r) 



W E |f ( k ) T E e -^G« (k, zfc m ) 2 



1 



g/(gO 



so that 



(C50) 



J P dr y £ i e- ik ^e i, ^ T fo(T j ,T)V(T j ,T). (C51) 



- 4 E f£ « T E G(2) ( k < ik ^ G(2) ( k > 

k ik m 



Using our previous definition, Eq.(C31), for f n {k2), the 
third term in Eq.(85) can be rewritten as 



U (T 



', ( v j E ^rQLi)GW(k 1 )GW(k 1 +i qn )GW(k 1 + qi +iq n ) 



x 3C/ 



1 



1 



J7d 



1 



SP 1 - %Xo(9i) 1 - %Xo(<Zi + ' 1 + ^Xo(«i) 1 + ^Xo(<Zi + »«„) 



x ^E/"( fc2 ) \G {1) (k 2 + qi+iqn) + G^(k 2 - qi )\ . (C52) 
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The sum over fc 2 is the sum of two convolutions and can be written as 

§ £ Uk 2 ) [G^(k 2 + q 1+ iq n ) + G {1 \k 2 - gi )] = ]T ( e *(fli+w»)-I + e ~^ 1 ) / n (l)G«(-T) 

k 2 1 

= ^Jo dT {^ T ^~ i(qm+qn)T + e-^e^) /„(!•„ r)G«(-r,, -r) (C53) 
= J* dr (-e-^+i^r + e i qm r^ J- e-^MTi^G^iri, -r) , 

where we have used f n (-Tj,T) = -/ n (rj,r) and G (1) (-r 3 , -r) = G (1) (r,-, -r). Using /„(r,r) = #„(r)/i(r, r), 
Eq.(C38), we have 



T 

TV 



E/«( fc 2) [G ( %2 + <Zi+i3„) + G (1) (A;2-<Zi)] = 

/ dr (. e ^ + e^) ftl (r)j;e- i *' r 'l 1 (r j) T)G( 1 )(r j) -T) (C5 4) 
= — f dr (- e -^ T + e- j(9 '" +9 " )T + e 4(9 '" +9 " )r -e^ T )^e^ qi - r ^(r„T)G (1) (r„-r). 

] 

assuming g„ ^ 0. Inserting this expression into Eq.(C52), we get 

™ fcl,9l 21 

X ( 3C/sp l- % X o(9i) l-%Xo(5i+«5n) + ^1 + ^X0(91) l + ^Xo(«i+*«n)) 

x I dr (- e -^ r + e- 4 ^ + ^) T + e J ( 9 ™+9") T -e J «'" r )Ee- 4C ' 1 - r ^(r J ,T)G (1) (r J ,-T). (C55) 

Using the definitions 

G«(fc! + qi ) = G^(h iq n ) (C56) 

and 

Hn[qi) = { WSP I - 4xo( qi ) l-%*o( 9 i+«fn) + ^l + ^Xo( ? i) l + ^Xo( 9 i+«Z«)) 

x f dr (e^m+q n )r + e -*(q m + qn )T _ ^q m r _ e -iimr^ ^ e~ icl1 ^ h(rj , t)G^ (rj , -t) , (C57) 

J 3 

Eq.(C55) reads 

fcl 1 



where the sum over was written as a Fourier trans- a function of k\ or fci + i(/„, or a sum of either, because 
form. Unfortunately, since this transform does not give the dependence on iq n in H n { q i) does not factor out as 
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a sum of exponentials, the sum over k\ does not have 
the form of a convolution or a sum of convolutions and 
has to be done in the form Eq.(C58) for each different 
frequency iq n . 

Now, some work remains to be done before one can 
do the Fourier transform over 1 in Eq.(C58). First, 



Gn \r,r) is explicitly given by 



1 



ik m + iq n 

(C59) 



so that changing summation variable above or using 
Eq.(C23), we obtain 

(r, r) = e iq " T G^ (r, r) . (C60) 

Eq.(C58) thus becomes 

E d ^)G {2) (^k m )G^{^k m + iqn ) 

iq n z iv Ohj X 

X E/ dTe lk V' fc+? " )T G (1) (r,-T)ff n (r,T). 

(C61) 

There is however a problem if we directly use H n (r, t) 
in Eq.(C61). From the definition (C57), one notices that 
H n (q,iq m ) is peaked at the frequency q. m — ~q n because 
of the factor 

(C62) 



1 - % 1 xo(q,«'?m + iq n ) 



in the spin part. Therefore, H n (r,r) has oscillations in 
r at the frequency q n and it is necessary to refine the 
grid in r when q n increases. H n (q 7 iq m ) is also peaked at 
q m = because of the factor 



f 



r xo(q,«5m) 



(C63) 



It is therefore not sufficient to express H n (r,r) using a 
function shifted in frequency multiplied by an oscillating 
function. By doing that we would reduce at most the os- 
cillation frequency by half. Instead, we express H n (r,r) 
as a sum of two functions, one peaked at q m = and the 
other, at q m = —q n . Then we will be able to apply a 
translation in frequency to the latter and factor out the 
oscillating part. First, we write 

1 1 

i - ^xoCq, iq m + iq n ) l - %xo(q, iq m ) 

i xo(q^gm) 

xo(q, iq m ) - xo(q, iq m + iq n ) 1 - ^x (q, iq m ) 
l 

Xo(q, iq m + iq n ) - xa(<\,iq m ) l 



XcrygMgm + iqn) 
■ %xo(q,i?m + iq n ) 



(C64) 



If we define 



D(q,iq m ) = 



J dr (e^ + e- i ^)J2e- i ^h(r j ,r)G^(r j ,-r) . 



the spin part of H n (q, iq m ) reads 
D(q, iq m ) - D(q, iq rn + iq n ) 



+ 



xo (q, iq m ) - xo (q, iq m + iq 

D(q, iq m + iq n ) - D(q, 

1 Qm ) 



Xo(q, iq m + iq n ) -xo(q,«9m) 



(C65) 

x S p(q, iq m ) 
Xs P (q, iqm + iqn)- 

(C66) 



Note that, since D(q, iq m ) and xo(q> iqm) are even func- 
tions, the function 



D(q, iq m ) - D(q, iq m + iq n ) 

xo(q, iq m ) - xo(q, iq m + iq n ) 



(C67) 



is undetermined when q m = —q n /2, which happens when 
n is even. However, we know that H n (q, — ^) vanishes 
so that one can assume an arbitrary value for Eq.(C67) at 
that point. Numerically, it is better for that factor to be 
smooth. This is achieved by using a simple interpolation 
to fix that value. 

Now, if we define the function 



I n {<l,iqn 



3t/„. 



D(q,iq m ) - D(q,iq m + iq n ) Xo(q,^m) 



xo(q, iq m ) - xo(q, iq m + iq n ) 1 



L xo(q, iq m ) 



+U C . 



D(q,iq m ) - D(q,iq m + iq n ) Xo(q^gm) 

xo(q,*9m) - xo(q,«<?m + iq n ) 1 + ^xoCq,^™) ' 



(C68) 

which is peaked at q m = 0, we make explicit that 
H n (q,iq m ) is peaked at two values of frequency 



H n (q,iq n 



-7 n (q, iq m ) + I n (q, -iq m - iq n ) . (C69) 



Using the bosonic Matsubara frequency representation 
for /„ (q, t) we obtain 

H n (q,r) =Tj2e- lqmT H n (q,iq m ) 

iq m (C70) 
= -I n (q,T)+e^ T I n (q,- T ). 

To compute /„(q, r) from I n (q,iq m ) we use the same 
procedure as for the previous Fourier transforms over 
Matsubara frequencies. Using 



D inf (q,iq m ) = 

X0,m/(q,«9m) = 



d(q) 

(iq m ) 2 

c(q) 

(iqm) 2 



(C71) 
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the asymptotic form for J„(q, iq m ) is 



C / (q,"z™) = 3[/.. 



d(q) 



SP (i<?m) 2 - %c(q) 



,.„ H . (C72) 



Using a complex plane integration and the residue theo- 
rem, we get 



I™f( q ,T) = 3U s: 



d(q) 



n B (-z sp (q))e- z ^ 



2z sp (q) 
-n B (^ P (q))e z - (q)T 



'2iz c h(q) 



-n B (iz ch (q))e 



»«ch(q)T 



(C73) 



where ns{z) is the Bose- Einstein distribution and 



>(q) 



fc(q). 



(C74) 



or, using ub(—z) 



C / (q,r) = 



n B (z), 



- 3U sp ^^n B (z sp (cO)e^^ cosh 
-^^^(i^lq))^^^' cos 



t- 2 sp (q) 



2 _r ) z ^(q) 



(C75) 



Using that result we write, as before, the transform in 
such a way that it converges quickly 

I n (q,r) =T£V i *» T [J„(q,*g m )-4 n/ (q»i9m)] 

+ C'(q,r). (C76) 



"In 



There is however one last difficulty. The asymptotic 
form Eq.(C72) approaches Eq.(C68) only when both q m 
and q m + q n are large with respect to the bandwidth. 
Thus, the frequency range of the finite sum in Eq.(C76) 
has to be chosen such that this condition is satisfied. If 
q n is positive, the range of negative q m must therefore 
be extended to make sure that, at the cutoff, q m + q n is 
large. 



Coming back to expression (C58), using the results 
(C60) and (C70), the sum over 1 becomes 

1 

J dT^ r e- l{k ™ +q " )T G {1) {v,-T) 



x [-J„( r , T ) + e ^ r /„(r,-r) 



E / dT 



■ e ik-r e -ik m T 



-e-*«»' G (i)(r,-T)I n (r,T) 

+G^(r,-r)I n (r,-T) . 

(C77) 



For q n = 0, expression (C52) vanishes, as explained at 
the end of appendix B. 

To reach reasonably low temperatures without seeing 
any finite size effect, we use a system of 512 x 512 sites 
and 8192 frequencies. If we were using a brute force 
approach to compute Eq.(C52) we would have to sum 
(8192(512 2 )) 3 = 9.9 x 10 27 terms. Assuming we could 
sum 10 9 terms per second, it would take about 300 billion 
years to calculate Xj*j (*9n) f° r onc value of q n and thus 
about 30 000 billion years for a hundred values of q n . 
Using the approach described above, this calculation is 
done in less than 2 days. 



Appendix D: Choice of Matsubara frequencies 

While all external frequencies are obtained at the 
same time from the last FFT for the bubble and Maki- 
Thompson like term in Xj x jx (*9«)> this is not the case for 
the Aslamasov-Larkin like term, Eq.(C52), even in the 
form of Eq.(C58) that makes maximal use of fast-Fourier 
transforms. Hence, this term cannot be calculated for 
thousands of values of q n , or for the same number as 
that used in internal Matsubara frequency sums. There- 
fore we have to compute Eq.(C58) for a reasonable num- 
ber of carefully chosen frequencies. Assuming that most 
of the information is in the low frequencies and that, 
as their magnitude increases, it becomes less important 
to include all the high frequencies, we use the following 
non-uniform Matsubara frequency index grid. That grid 
consists of subintervals within which the Matsubara fre- 
quencies are equally spaced, with larger space in between 
Matsubara frequencies at large frequency. The spacing 
between frequencies in different subintervals increases by 
powers of 2. 

First we define 

N ~ 2 r : the last frequency index (cutoff), taken as a 
power of 2 (r integer), 

m: integer that determines how dense the grid is. A 
large m gives a low density. < m < r, 

Ni — |sr : Ni + 1 is the number of adjacent frequencies 
close to n = 0, 
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N 2 — '■ number of frequencies in each subinterval 
with a fixed spacing between Matsubara frequencies, 

N = Ni + mN 2 + 1 : total number of frequencies. 

Then the indices of Matsubara frequencies on the grid 
are given by the following algorithm 



0,...,^!-!, 



Aq + 2 l i+ l mod{j - N U N 2 ) + iVi(2'i - 1) , 



lj = floor , j = N 1 ,...,N-l. 

(Dl) 

For example, taking N — 256, m — 5, so that N\ = 8, 
N 2 = 4 and a total number of frequencies N = 29, we 
obtain the following indices: n = 0, 1, 2, 3, 4, 5, 6, 7, 8, 
10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 64, 80, 96, 112, 
128, 160, 192, 224, 256. 

With this kind of grid we greatly reduce the number 
of frequencies for which we have to calculate Xjlj x with 
Eq.(C58) while retaining the essential information. Note 
that this also speeds up the analytical continuation with 
our maximum entropy method described in appendix F. 



Appendix E: Fourier transform of a cubic spline 



Assume we have the following integral to do 



/(*) 



dx g{x)e 



ikx 



(El) 



but that we know only N + 1 discrete values of g(xi). Let 
us approximate g(x) in the interval using a cubic spline 
S(x) defined as 



S(x) 



S\(x) x < x < x\ 
S 2 (x) x\ < x < x 2 

Sn{x) xn-i < x < xjm . 



(E2) 



where the S n {x) are cubic polynomials, with the condi- 
tions 



S n (x n -i) = g(x n -i) 
S n { x n) = g(x n ) 

S'n (.Xn- l)=S' n _ 1 (X n - 1 ) Tl > 1 

S%{x n -i) = S^fa-i) n>\ 
S[(x Q ) = g'(x ) 
S' n (xn) = g'{x N ) , 



(E3) 



defining the AN equations necessary to determine the AN 
coefficients of the spline. The integral (El) becomes 

N rx n 

dxS n {x)e- ikx . (E4) 



Integrating by parts we obtain, 



n=l 



-ikx 



i 

ik 



N 



E 

71=1 

e -ikx oSl ( Xo 



S n (x) 
dx S' n (x)e- lkx 

-1 

e- ikx " S N (x N ) 



(E5) 



ik 



+ 



1 N 

-T 

ik ^ 

n—l 



dx S' n (x)e 



-ikx 



where we have used the continuity of the spline at the 
points x n to eliminate all the intermediate terms in the 
first sum. Now, if we integrate by parts in the second 
term and use the continuity of the derivatives S' n (x) at 
the points x — x n , we obtain 



f(k) 



*kxo Sl ( 



Xo) 



S N {x N ) 



ik 



+ 



e- lkx °S[{x ) - e~ ikxN S' N {x N ) 



{ikf 



N 

V 



dxS'^{x)e 



-ikx 



(E6) 



Doing it one last time, using the fact that we also have 
the continuity of the second derivatives S'^ (x) at x n , we 
finally obtain 



/(*) = 



e' lkxa Si(x {) ) -e~ ikxN S N {x N ) 



+ 



+ 



ik 

- ikxo S[(x ) 



-ikx 



N S' N (x N ) 



(ik) 2 
>S>{{x Q )-e- ik > 



'S%{x N ) 



(ik) 3 

-ikAx N - 



+ 



ik) 4 



y, ,s 'i 3 +i ( 



tkx n 



n=0 



(E7) 



where Ax = x n+ i—x n . The remaining discrete transform 
can be done using an FFT. The result of this transform 
will be periodic since it is discrete, but because of the 
factor 1/ (ik) 4 in front, it will only have a relevant contri- 
bution for small k. This periodicity in the sum therefore 
produces only a very small noise at high k, the remain- 
ing discretization noise, that decreases with increasing k. 
Note that expression (E7) is valid only for k ^ 0. For 
k = 0, the result is simply the sum, over all subintervals, 
of integrals of cubic polynomials with their respective co- 
efficients. 

Finally, note that we have chosen to fix the value of the 
derivatives at the boundaries to complete the system of 
equations (E3) defining the spline. But other interesting 
choices are also possible. For example, if we know the 
coefficients of the terms in l/(k 2 ) and l/(fc 3 ), i.e. the 
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numerator in the second and third terms of Eq.(E7), fix- 
ing those coefficients is a good alternative. In the case 
where k is a frequency, this choice is convenient because 
we often know the high frequency expansion from sum 
rules. 



Appendix F: Analytical continuation for the 
conductivity 

Let us start by rewriting the spectral representation 
for the Matsubara current- current correlation function 



Xj 



.j.(*«n) = y ^ 



du x'LM 



Jxjx 

u - iq n 
du x'LiA") 



. [du [du 
J ir u 2 + ql J IT 



duj wxlixM 



u 2 + ql 



Since x'j x j x ( w ) is °dd, 



Xjx 



w 2 + <? 2 

The real part of the optical conductivity is 



<7» 



so that we obtain 



Jo 



7T w 2 + q- 



-a r (u) 



(Fl) 



(F2) 



(F3) 



(F4) 



using the fact that the integrand is even. This is our 
starting point. The objective is to obtain the real part of 
the conductivity, on the right-hand side, from the Mat- 
subara expression for the susceptibility on the left-hand 
side. Most analytical continuation is done for imaginary- 
time data, but not in our case. 75 

Suppose we could be satisfied with cr r (u) on a discrete 
set of points Uj. We can use a numerical integration 
method to approximate the integral (F4), which would 
then have the form 



Xj.j.(*«n) ;i: E K 



(F5) 



where crj = <J r {uj) and K n j is a N qn x JV^ matrix, N Qn 
being the size of the vector Xj x j x (iQn) an d N u , the size 
of the vector crj . Now, crj is the quantity we want to de- 
termine. If N qn — N u then <rj is completely determined 
by the linear system (F5). However, the matrix K n j is 
ill-conditioned so that a small noise in Xjxjxi^ln) would 
result in a very noisy solution crj. Also, N Qn is gener- 
ally smaller than N u , the number of real frequencies for 
which we want to determine <J r (u). We therefore need to 
include more information in the problem to find a unique 



<rj. The way to do this is to use a maximum entropy ap- 
proach. In this approach, we minimize the function 



where 



* 2 = E 



aS , 



Xj*j x ( i >ln)-T l j K njC7 r j 



(F6) 



E 



(F7) 



measures the deviation of Xj^jxi^Qn) = J2 3 Knj&j with 
respect to Xj x j x {iQn), e n being an estimate of the error 
of Xj x j x (iQn) with respect to the "exact" function. S is 
the differential entropy defined as 



Jo 



du> <j r (to) In 



m(w) 



(F8) 



N Qn , so that 



where m(u>) is called the default model. The value of 
a can be chosen according to different criteria. As is 
often done, we choose it such that x 2 
\Xj x j x {iQn) ~ Xj x j x {iq n )\ be equal to e„ on average 

Errors in the numerical evaluation of the integral (F4) , 
i.e. in the definition of K n j in Eq.(F5), are equivalent to 
having larger errors in the data Xj x j x (^Qn)- Therefore, 
when Eq.(F6) is minimized to find a solution erj, this 
could lead to large errors in crj with respect to the op- 
timal solution because the inversion of expression (F4) 
is an ill-conditioned problem. It is thus clear that the 
error we make by replacing Eq.(F4) by Eq.(F5) must be 
smaller than the estimated error e n on the original data 
Xj x j x {iQn)- That is why we need a very accurate numer- 
ical integration technique to define K n j. Because we use 
the spectral representation in Matsubara frequencies, the 
weight function oj 2 /(u> 2 +(?„) in the integrand of Eq.(F4) 
is simple and can be integrated analytically. Hence, if we 
use, for example, a polynomial approximation for <t t (uj) 
in a given interval [uj_i,ujj], the integral can also be 
done analytically in the interval. If we use a good piece- 
wise polynomial approximation for <j t (uj), then we can 
evaluate the integral (F4) precisely. 

Maybe the most efficient approach to integrate Eq. (F4) 
in one dimension is Gaussian quadratures. However, in 
the latter, both the weights and the grid points depend 
on the weight function, which in our case depends on 
q n . We would therefore need a different grid in u> for 
each frequency iq n . This is not possible because we can 
only search for a unique vector crj defined on a unique 
grid uij. Another very efficient way of doing the inte- 
grals (F4) is to model a r (u) using a cubic spline. This 
approach allows to use a fixed grid and is very precise be- 
cause cubic splines are very good to approximate smooth 
function, which is the case of <J r (u). It seems therefore 
the best approach for our problem. However, to be able 
to perform the integral (F4) over the whole frequency 
range [0, oo, and at the same time reduce the number of 
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frequencies in the grid, which helps speed up the min- 
imization of Eq.(F6), the spline we use is divided into 
two parts, a low frequency part that is a cubic spline in 
oo and a high frequency part cubic in u — 1/uj. Then, to 
make the spline linear system well-conditioned, we use a 
grid that is uniform in oo in the low frequency part and 
uniform in u in the hight frequency part. Finally, inte- 
grating analytically in a piecewise manner, keeping the 
weight function oo 2 /(oo 2 + q^) intact in the integrand is a 
great advantage as the temperature decreases since this 
function then becomes sharper and sharper, and is thus 
increasingly difficult to integrate numerically. Of course 
the conductivity <J r (oo) itself becomes also sharper as T 
decreases and we have to adjust the grid to resolve its 
structure. But the numerical integration is still much 
easier and precise with this approach. We describe be- 
low how the K n j is defined and also our choice of grid in 

00. 

We start with the following representation for cr r (oo), 
a r (oo) = 

!Sj{uo) , OOj-! < 00 < LOj , 1 < j < N 

8j(±), uoj-i < oo < ooj , N < j < N + M , 

(F9) 

where Sj(x) — ajx 3 + bjx 2 + cjx + dj, with the conditions 
s jK-i) = . 

s 'j(uj-i) = Sj-iK-i) , (F10) 

s "K-i) = s 'j-i K-i) ' 
«i(0)=0, 

for j < N, while 



sn+i [ — ] = CF N . 
, 00n , 



SN+1 



00n+1 

1 



'JV+1 



oo N 

1 

OON 



and 



OOj-! 



= "N+l i 

= s' N (oo N ) 
= s' N ( w Jv) 

r 

J 3-l ' 



(Fll) 



OOj-! 



S J-1 



OOj-! 



OOj-! 



(F12) 



ds N+M (£) 



for N + 1 < j < N + M. All the derivatives noted with 
' and " are taken with respect to oo. The second condi- 
tion for j = N + M and the last condition in Eq.(F12) 
make sure that there is no constant and no 1 /oo terms in 
sn+m(1/oo) so that the integral of the spline converges. 
Note that this spline is not physically correct for oo — >■ oo 
since the moments of <J r (oo) are not defined. But this does 
not have any significant importance in the numerical so- 
lution if oon+m-i is chosen large enough. The advantage 
of this spline is that it is very easy to implement. 

Let us now define the kernel K n j in Eq. (F5). For the 
low- frequency part of Eq. (F4) we have 



["» doo oo 2 



oo z + qi 



9 N 

= -E 



dw - .2 + n 2 "3 

-l oo + q n 



Sj {oo) , 



N 



T j=1 J Ui -l 



, cijoo 5 + bjoo 4 + Cjoo 3 + djoo 2 
doo — 



oo 2 + q 2 



(F13) 



where oo — 0. For q n ^ 0, we obtain 
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Wi- 1 



+ 



oo 2 q 2 



y^(ql + oo 2 ) 



+ 



and, for q n = 0, 



oo — q n arctan 



oo 



(F14) 



(o) = I E / J doj ( a ^ 3 + b ^ 2 + c ^ + d i) 

9 N 



o N / 4 , ,4 , ,3 , ,3 

2 ^ / OOj - OOj-! „ "j - Uj-! 

4 aj+ 3 



a; 2 — oo 2 _ 1 

+ 3 o J C 3 + ("j - Wj-l) rfj 



The high frequency part of Eq. (F4) is 



Using 



w = 



doo = -du , 



(F15) 



(F16) 



(F17) 
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it becomes 



u M 



du 



1 



ir J u z + u 4 q2 

where u M — 1/uj n . Now, with Eq.(F9), we have 

M 



U 



(F18) 
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(F19) 



where 1 < j < M, uq = 0, and we have used a different 
notation for the coefficients to match their indices with 
those of grid points in u, which are indexed in order of 
decreasing w. For q n ^ 0, we obtain 



M 



and, for 
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(F20) 
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(F21) 



Note that 71 = and Si — 0, so that this expression has 
no problem with u = 0. 

The spline coefficients for the complete spline are ob- 
tained at the same time for both the low frequency and 
high frequency parts using the conditions (F10), (Fll) 
and (F12). Because the spline has two parts made of 
polynomials in w and u = 1/lu, one must define the spline 
linear system for that particular case, but this is a rather 
straightforward task. Note that, for the high frequency 
part, the condition of continuity of the first derivatives 
with respect to u or ui are the same, but that is not true 
for the continuity of the second derivatives. However, we 
do not observe any loss of accuracy if the continuity of 
the second derivatives with respect to u instead of oj is 
used. 



The conditions (F10), (Fll) and (F12), put in matrix 
form are written as Av — a, where v is a vector contain- 
ing the spline coefficients and a is a vector containing 
the values crj (that are repeated) and zeros otherwise. 
We want to obtain the coefficients as linear forms of the 
vector a r . For that we have to invert the matrix A and 
first obtain a linear form for v. The first two lines of 
the condition (F10) tell us that a contains the elements 
crj repeated twice each, except for erg, and the other el- 
ements of a are zeros. So if we define the matrix Pao- 
that sum up the pairs of columns in A^ 1 corresponding 
to the same crj and removes all columns that correspond 
to zeros in a we get 



v = A- x P Aa o r 
= T va a r 



(F22) 



Then the vectors formed with the coefficients are given 
by expressions like 



P a v, 



P h v. 



(F23) 



where P a extracts a column formed of all the coefficients 
of the cubic terms, P b the coefficients of the quadratic 
terms, etc. Now, expressions (F14) and (F15) have the 
matrix form 



*f.i. fa") = " ( K > + K n b + K nC + Kd\, (F24) 



which becomes, using Eq.(F23) and Eq.(F22), 



A. = ~ [KnPa + ^£ A + K°P C + R d n P d ^ T va O r . 

(F25) 

Similarly, if we define the projectors for the vectors a, 
/3, 7 and S, then Eq.(F20) and Eq.(F21) take the form 



X 



hf 



(iq n ) = - ^P a +KHP fs +KZP ( + K S n P S jT v<T a r . 

(F26) 



Summing Eq.(F25) and Eq.(F26) we obtain 

Xj*jMn) = K n a T , (F27) 

with 



K„ 



K a n P a + K b n P b + K c n P c + K d n P d 



+ K«P a + K?P fj + K^Pj + K s n P^j T va . 

(F28) 

To conclude this section, a few other points must be 
addressed. First, the form of expression (F14) becomes 
numerically unstable when q n /uj are large. For exam- 
ple, the first two terms in the high q n /uj expansion of 
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in K°: cancel out the terms 



■\uj 2 and 

When q n /ui increases, the magnitude of those terms 
becomes much larger than K° itself so that if one com- 
putes numerically the three terms of separately and 
then adds them, the finite precision error becomes larger 
than the true result. To overcome this problem one sim- 
ply has to use the large q n /u> expansion of K% starting 
at a certain cutoff, instead of directly the form appear- 
ing in Eq.(F14). The same kind of cancelation appears 
in the other terms of Eq.(F14) so that the large q n /u> 
expansions must be used for those terms as well. In the 
case of Eq.(F20) it is when q n u is small that some sim- 
plifications occur when replacing the expressions by their 
expansions, which will also improve the accuracy of the 
numerical result. 

Second, since we invert the spline matrix A to define 
T VIJ in Eq.(F22), it is preferable that this matrix be well- 
conditioned. For that purpose we define our grid to be 
uniform in uj for the low frequency part of the spline and 
uniform in u for the high frequency part. The grid in u 
is 



1 



1 



u = 0, 

If Uj — Uj-i is constant we have 



W/V+M-2 — 



Wjv+M-1 



Wat+M-3 — 



UJ N 
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UJ N 
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(F29) 
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(F30) 



or 



Muj n 
M-l ' 



— 



Muj n 
M-2 ' 

un+m-i = Mujn ■ (F31) 



Also, to have a density of points ojj that varies con- 
tinuously when we change from high to low frequency, 
we assume that lon+i — ^>n = — wjv-i- Defining 
Auif — ojn — wjv-i, we have 



Muj N a 
w - J -uj N = Auj lf , 



so that 



and, if i 



M 



Awif 



+ 1 



-i is constant for j < N, 



Auji f = 



N 



(F32) 



(F33) 



(F34) 



so that M = N + 1. Note that when choosing uj n and 
N, we have to check that the last frequency ujn+m-i = 
(N+1)con is large enough while not so large that it would 
uselessly make the calculation heavy. To further improve 
the conditioning of the matrix A, we use normalized fre- 
quencies lo'a — ujj/ujN so that uj' n — 1 and u' M = 1. In 



the j th interval, if a'-, c' et d' are the coefficients in 
the normalized grid, we have 



bjUJj 



CjUJj 



= a. 



+ c'j 



LO.j 



J uj n ■'ujn 



+ d' 



(F35) 



and therefore, 



3 

4, 



UJ N 



, dj = d\ . (F36) 



Finally, for the high frequency part of the spline, we have 



aj = uj Naj , 



5,=^. 
(F37) 

To end this section, we comment on the differential 
entropy and the minimization procedure. For the differ- 
ential entropy (F8), we use a default model that is almost 
fiat in the region where <r r (uj) is expected to have its main 
structure, and that decreases gradually to very small val- 
ues for frequencies much larger than the bandwidth. This 
ensures that the solution that we find is as unbiased as 
possible. The minimization of Eq.(F6) is performed us- 
ing a Matlab routine called fmincon, which uses a Trust- 
Rcgion-Reflective algorithm that has been proven quite 
efficient with not as much tendency to get trapped into 
local minima as other optimization routines. Our proce- 
dure is to start with a very large value of a such that the 
minimization process gives a solution very close to the 
default model m(uj) (the minimization of Eq.(F8) alone 
has in fact a solution proportional to m(w)). Then, a is 
decreased and a new optimal solution is found, using the 
previous solution as a starting point in the optimization 
routine. This step is then repeated until \ 2 ~ N qn or x 2 
does not decrease anymore when a is reduced. Using an 
augmented lagrangian method, we also include inequal- 
ity constraints to restrict the roughness of the solution 
<j r (ujj). This roughness appears at some point in the pro- 
cedure when we try to make K n a r closer to Xj^j^i^Qn)- 
It is related to oscillations present in K n j as a function 
of the frequency index j for a given n. Those oscillations 
are in fact the price to pay to work with an accurate nu- 
merical integration method since they are present in all 
approximations more sophisticated than a piecewise lin- 
ear function for a T (uj) in the numerical integration (think 
about the unequal weights in a Simpson 1/3 or in gaus- 
sian quadratures). In fact, the oscillations in our K n j 
have an extremely small relative amplitude, but they ap- 
pear greatly amplified in the solution <r r {ujj) when the 
relative distance IXjWx — K n cr r \ becomes very small. 
The link between oscillations in <r r (ujj) and oscillations 
in K n j for n fixed is clear because they are correlated in 
the two functions. So they are not related to noise in the 
data XjxjAiln)- However, those inequality constraints 
are not absolutely necessary to obtain good, quantita- 
tive, results with our approach. They just ensure that 
a r (ujj) is a smooth function. 
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